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ABSTRACT 

Two related models are proposed to explain aspects of 
the connection between seismicity , water depth, loading 
history, and hydrologic conditions at an artificial lake. 
Both models involve fluid loads on layered half spaces. The 
response in each case is most rapidly determined with 
Fourier transform techniques. The models are laterally 
homogeneous, but in spite of this simplification predict a 
gratifying number of observed seismic features. A model of a 
load on an elastic half-space predicts many features, but 
does not provide time dependent predictions. A similar 
calculation for a permeable layered structure predicts that 
initial seismicity will occur directly beneath a reservoir 
during the first filling. However, the region most subject 
to seismicity changes in pcsition and time as the lake depth 
varies and may be offset from the taney The response of a 
permeable structure can be expressed as a product of two 
terms; one evaluated at the particular time of interest and 
the other a convolution integral containing the rate of 
change of depth. 

Examples of the calculations for various loading 
histories are given in two and three dimensions. For a two 
dimensional lake model and an assumed filling history, it is 
shown that seismic risk is a function of time and position 
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with respect to the lake. Maintaining a constant water depth 
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allows the region to stabilize slightly by permitting lobes 
of high pore pressure to diffuse away. It is also shown that 
a partial emptying and refilling sometime after the initial 
filling can cause extremely large stresses in the rock. In 
some instances it is felt that draining and refilling may 


cause an increase in the regional seismicity. 
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The filling of large artificial lakes has often been 
associated with an an apparent increase in the seismicity of 
the region. This work will investigate the seismicity fron 
the theoretical aspect of examining stresses produced in the 
rock by the application of a surface load - the lake. A 
mathematical approach has many difficulties, since the 
application of simplified theoretical models to the 
complicated Earth rarely leads to results which are 


conclusive in every case. 


The literature on seismicity around lakes is extensive 
with the current activity directed toward the aim of being 
able to predict whether the filling of a certain lake will 
increase the region's seismicity. This is called the 


“predictive problem". 


Perhaps the most important consequence of lake 
seismicity can be expressed in terms of the life and 
property damage which may occur near epicentres of 
magnitudes greater than 6 on the Richter scale, several of 
which have been detected close to large reservoirs. The 


descriptions of these earthquakes can be found in several 
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excellent review articles (See Simpson (1976), Gupta and 
Rastogi (1976)) and need no detailed documentation here. At 
Koyna (see Figure 1.1) 200 lives were lost and a large 
amount of property damage was caused by a magnitude 6.5 
earthquake within 10 km of the lake. It is generally 
conceded that dam engineers can design structures to 
withstand extremely large earthquakes safely so the chances 
of a dam failure are minimal. They would however like to 
keep the cost down by knowing whether an earthquake is 
likely and, if so, its probable size. An assessment of the 
risk to the surrounding structures and population is also 


part of the predictive problem but is not discussed here. 


Simpson (1976) has presented an excellent summary of 
the seismicity around reservoirs. What is obvious from his 
list of ‘large* reservoirs is that not all of them are 
associated with earthquakes. The definition of a large 
reservoir given by the Ccmmission On Large Dams is one that 
is deeper than 100 metres and contains a _ volume of 


approximately 1.2 billion cubic metres. 


This thesis will concentrate on some of these ‘large’ 
dams since the weight and pressure head at the at the bottom 
of these is the greatest. Large dams also seem to be _ those 


most subject to noticeable seismic changes. 


In the World Registry of Dams (I1.C.0.L.D. 1973) there 


appear 275 dams with heights larger than 100 metres with a 
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Figure 1.1: 


A world map showing the positions of most 
of the lakes with significant seismic 
changes after filling. Rangely and Denver 
are sites where tests on the dependence of 
the seismicity on fluid pressure were 
conducted. 
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further 149 planned or kteing constructed. Of these 424, 93 
are over 150 m@ and 38 have volumes in excess of 10 billion 
cubic metres. It is apparent that the size of dams, as well 
as their number, is rapidly increasing. Ihe tallest is Nurek 
at 317 m high. Bratsk contains a volume of 169 billion cubic 


metres of water. 


Table 1.1a (modified from Simpson (1976)) lists the 
largest dams and Table 1.1b the deepest. There appears to be 
ho consistency between either the water depth or volume and 
the seismicity. This indicates that neither depth nor volume 
are in themselves sufficient to change the seismic pattern. 
Rothe (1973) has indicated that depth appears more important 
than the total volume and that the activity is greatest in 
reservoirs deeper than 100 metres. This is not evident fron 
the tables 1.1a, 1.1b and the 100 m depth criteria now 


appears meaningless. 


In many cases lakes are situated in areas of relatively 
high tectonic activity. Even without the presence of the 
lakes, a large number of events can be expected. The problem 
of trying tc separate parameters which indicate a change in 
the seismic level in the region before and after the filling 
is poorly treated in the literature. One problem to be 
considered is how far a seismic region extends. Gough and 


Gough (1970b) examined the seismic history for the area 
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TAB hemlars HIGHEST DAMS 
DAM NAME LOCAETON DATE SEIS- TYPE DAM RES. 
MIC HGT VOI. 
m m°x10° 

NUREK USSR Te * ER 35157 O00 
GRANDE DIXENCE SWITZERLAND 62 PG 285 400 
INGURI UeShoie: G VA Dy ia2 1100 
ROSSELLA 4 se 65 TE 2015 Le) 
VAJONT etal, 61 * HA Il Bho) S: 
MICA CANADA 74 TE Dea? 2086s) 
SAYANSK US SPs G VA Zhe 30300 
MAUVOISIN SWiLt Ze RLAND 57 VA Dey" | 180 
OROVILLE USAre CALLE 68 * ee, 236 14298 
CH URKEY, USSR C VA 255 2780 
ESMERALDA (CHIVOR) COLOMBIA G75 ER 230 815 
BORUCA COSTA-RICA P80 ER 229 617 :0:0 
BHAKRA INDIA 63 PG 226 9868 
HOOVER (L. MEAD) DSA, AREZ.NEVeD 6 VA 224. 36703 
MRATINJE YUGOSLAVIA C VA 220 880 
CONTRA SWITZERLAND 65 * VA 220 86 
PALLA COLOMBIA iP ER 220 11000 
DWORSHAK DAM USA, IDAHO i383 PG 2A9 42:76 
GLEN CANYON (POWELL)USA, ARIZ. 64 VA 246 33305 
TOKTOGUL USSR C VA 205 L9500 
DANIEL JOHNSON CANADA 68 MV 204 Wa8od: 
AUBURN USiAn 9 sGAdsl 0 C 1g 210 283F 
LUZZON#: SWLIZEREAND » 63 VA 208 87 
KEBAN TURKEY 74 ER 207 31000 
HIGH MOUNTAIN SHEER USA, IDAHO P VA 204 4 
MOHAMAC REZA CHAH IRAN 63 VA 2038 8340 
URANCANCHA PERU e VA 203 70s. 
ALMENDRA SPAIN 70 VA 202 2649 
ROSS HIGH USA, WASH. 2 VA 201 42638 
GRAND-MAISON FRANCE ie VA 200 205 
RAZA CHAH KABIR IRAN 13 VA 200 2900 
KONDJILA YUGOSLAVIA P VA 200 416 
NE eS WNC 1d TE =PEARTHELILL DATE - YEAR COMPLETED 

ER - ROCKFILL C - UNDER CONSTRUCTION 

PG - GRAVITY Ie - PROPOSED 

CB — BUTTRESS 

VA - ARCH 


- MULTI-ARCH 


TABLE 1b: 


DAM NAME 


BRATSK 


SAAD-EL-AALI (ASWAN) 


KARIBA 


AKOSOMBO MAIN DAM 


DANIEL JOHNSON 
KRASNOYARSK 
W.A.C. BENNETT 
ZEYA 

SANMEN HSIA 
CABORA BASSA 
USile ear Ml 
TANKIANGKOW 
HOOVER (L. 
LOUILOU 


MEAD) 


GLEN CANYON (POWELL) 


SAYANSK 
KEBAN 

MICA 

KENNEY 
FURNAS 
TOKTOGAL 
GURI 
ILTUMBIARA 
TAMAHARA 
KOLYMA 
TARBELA 
DICKEY 

SAC SIMAO 
MOSUL 
BHUMIPHOL 
GRAND COULEE 
GORDON 

HS INFENGKIANG 
NAGARJUNA SAGAR 
PATIA 
SOUAPITI 
MANICOUGAN 3 
NUREK 


NOTATION - 
LY PoE 


LARGEST RESERVOIRS 


LOCATION DATE SEIS— TYPE DAM RES. 
MIC Bele VvOL. 
m m-°x10°6 
USSR 64 PG IWSe lo927 0 
EGYPT 70 ER 111 164000 
RHODESIA 59 * MV 128 160368 
GHANA 65 ER 141 148000 
CANADA 68 MV 214 141851 
USSR 67 PG te ares 3 00 
CANADA 67 TE 183 70309 
USSR GS PG 115 68400 
CHINA 62 PG LO Jenko5.005 
MOZAMBIQUE 74 VA 171 63000 
USSR C PG 105 5:93.00 
CHINA 62 PG 12.09 5.1000 
USA ARIZ.NEV 36 * VA 22 1 23.67.03 
ZIARE P PG 37 335.000 
USAS w ARIZ). 64 VA 2i1.65 433.3015 
USSR C VA 242 ~ 31300 
TURKEY 74 ER 207 # 31.000 
CANADA 74 DE 242 24670 
CANADA 52 ER LOGS 22203 
BRAZIL 62 ER 127 20860 
USSR le VA 215 19500 
VENEZUELA 68 PG 1'OiGamelasi0 0 
BRAZIL P79 ER/PG 106 17027 
JAPAN P ER 116 16300 
USSR C ER 130 14600 
PAKISTAN 7h) ce Diy RReeltes | glassy, 
USA, MAINE P TE 1065-" 13048 
BRAZIL P TEPC) 20a 12500 
IRAQ P ER si 1250.0 
THAILAND 64 MiAchP Ge Side «1 22-00 
USA, WASH 42 PG 168 11975 
AUSTRALIA 74 VA 140.6 1 7A8 
CHINA 59 * PG LO See 500 
INDIA G TE Wee oll 165 
COLOMBIA P ER 220 LOO 
GUINEA P TE 121 gis 
CANADA C75 * TE 108) 10423 
USSR He, x ER 37 60400 
—- EARTHFILL DATE - YEAR COMPLETED 
—- ROCKFILL C - UNDER CONSTRUCTION 
—- GRAVITY P - PROPOSED 
—- BUTTRESS 
- ARCH 


—- MULTI-ARCH 
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around Kariba but no events were found. However, when the 
seismicity for the whole of southern Africa is plotted, as 
shown in figure 1.2, two seismic areas become apparent on a 
feature passing through Kariba. One of these is near the 
reservoir with all the events occurring after the filling, 
and the other in the Okavango. In both areas there exist 


large quantities of surface water on a tectonic feature. 


At the other extreme, Lake Kremasta was reported by 
Comninakis et. al. (1968) to show good correlation between 
the filling of the iake and the local seismicity. However, 
it is difficult to distinguish the epicentral distribution 
pattern before and after the fiiling using the same N.O.A.A. 


data bank. This is illustrated in figure 1.3. 


If we examine the seismicity within an area subject to 
the same tectonic forces (earth forces causing large scale 
deformations) several observations can be made at lake 
sites. In many cases the epicentres are usually very close 
to the lake - within 25 km and often much closer. Local 
arrays often place many events directly under the lake. 
Several examples of the epicentral distribution of the 
larger (magnitude greater than 3) events at Koyna, Kariba, 
Nurek, Oroville and Hoover are shown in figures 1.4 to 1.8. 
The positions and magnitudes of these events were determined 
from teleseismic data and may differ from those determined 
by local networks. The difference in position was about 0.2 


degrees at Koyna. 
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Figure 1.2: 


Epicentres plotted for Southern Africa 
from the N.O.A.A. Summary tape. Notice the 
linear pattern of events through Kariba 
and the Okavango Swamp to the south west. 
These are twc areas where there is a large 
amount of ground water. The intense 
distribution of events near Pretoria is 
probably related to the Johannesburg 
mining activity. 
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Figure 1.3: 


Epicentres from the N.O.A.A. Summary tape 
for the Kremasta region in Greece with the 
coast line and drainage pattern. (a) shows 
the distribution up to and including 1965, 
and (b) shows the same data after filling 
in 1965. By comparing patterns, it is not 
cbvious that a major change in the 
seismicity has occurred since filling. It 
should be emphasised that before 1965 
fewer events were recorded and their 
magnitudes were generally unknown. 
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Figure 1.4; 


Epicentres determined near Koyna in India. 
All the events shown occurred after the 
completion of the reservoir. Uncertainties 
exist in the exact position of any event 
by up to one-tenth of a degree. 
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Figure 1.5: 


Epicentres from the N.O.A.A. Summary tape 
for the Karika region in Rhodesia. All the 
events plotted in this regicn occurred 
after the completion of the reservoir. 
Notice the tight clustering of events at 
the deep end of the lake. Many more 
smaller events were observed by local 
stations but the smallest events plotted 
have magnitudes between 3 and 4 Mb. 
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Figure 1.6: 


Positions and baghitudes of events 
recorded before (a), and after (b), 
impoundment at Nurek in the USSR. Several 
events are seen close to the lake both 
before and after filling and it is 
concluded that there is no significant 
change in tectonic activity. Local arrays 
have observed migration of microseismic 
events towards the lake. 
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Figure 1.7: 


Epicentres near Lake Oroville in 
California up to and including 1968 are 
shown in (a). This shows a seismic gap 
within a 50 km radius of the lake. Since 
the filling of the lake,(b), the pattern 
has changed little. fhe position of the 
August 1975 magnitude 5.9 Mb event is 
shown by the star in (b). 
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Figure 1.8: 


The earthquake epicentres at Hoover dam 
are shown since the filling in 1936. 
Before this, only two events of 
undetermined magnitude were observed in 
the far NE corner of the area. The small 
number of events before filling is 
probably only a consequence of the 
shortage of instrumentation at that time. 
A cluster of events is observed near the 
dam wall. 
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In only a few cases have depths been determined. Local 
observations and the teleseismic data all indicate that the 
hypocentres are shailcw. Gupta et al (1970) have plotted 
positions determined from the local array at Koyna and the 
important results are shown in figure 1.9. The uncertainties 


in the determination are not given. 


Figure 1.9 does serve to show that the bulk of the 
events cccur less than 10 km deep but may extend to 30 kn. 
They are limited laterally to within 10 km of the centre 
axis of the lake. Reports of noise (Simpson 1976) during the 


earthquakes is indicative cf shailow events. 


Migraticn of seismicity has also been observed at some 
lakes. Simpson (1976) and Soboleva and Mamadaliev (1976) 
have indicated the arrays of Nurek show a migration of 
events toward the lake. They propose a model in which the 
increased pressure head below the lake migrates away from 
the lake. The pressure changes are delayed behind the 
elastic deformation by the diffusion lag. The region of 
failure moves toward the lake in such a model. Inherent in 
this model is the assumption that the elastic load falls off 
more rapidly with distance than does the fluid pressure; 


this will be examined later. 


The focal mechanisms (Sykes(1967,1970), Bufe et. al. 
(1976), Gough (1976) and an excellent article by Langston 


(1976)) observed at different lakes are consistent with 
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Figure 1.9: 


The array study at Koyna yielded the 
positions shown for the aftershocks. The 
events are clustering at the deep end of 
the lake and are relatively shallow. Most 
of the events were above 10 km deep and 
lay within 25 km of the lake. 

(Modified from Guha et. al. (1974)) 
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neighbouring tectonic faulting. At Kariba, Kremasta and 
Oroville normal dip-slip faulting is observed, while at 
Koyna, Hsinfengkiang and Hoover the mechanism was strike- 
slip. Jacob et. al. (1976) have reported a slight change in 
activity at Tarbela which is in a thrust faulting regime. 
Nurek has also had many small events in a thrust regime. The 
magnitudes of the main shccks near reservoirs have been as 
high as 6.5 at Koyna, 6.3 at Kremasta, 6.1 at Hsinfengkiang 
and so on down. There would appear to be no upper limit to 
the magnitude since recent work suggests that the lake acts 
only as a trigger. The magnitude is then only limited by the 
stresses available in the region. The fact that there have 
been only a few events cf magnitude greater than 6 near 
jakes can be explained in two ways. Firstly, the lake has to 
ke placed in a suitably highly prestressed region; and 
secondly, the larger events are less frequent than the 
smalier ones. Scme authors feel that the magnitude of the 
Main shock is limited since the anomalous stresses due to 
the lake are localized and the rupture length of a large 
earthquakes is hundreds of kilometers. This argument ignores 
the possibility of a locally induced rupture extending into 
an area in which there exist sufficient stresses for the 


fracture to continue to propagate. 


The magnitude-frequency distribution of most 
earthquakes sequences follow the empirical rule 


Log N = A-bM 
10 
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where N is the number of shocks of magnitude greater than or 
equal to M. "A" and *b* are constants determined for the 
region. The constant A is uncertain since it depends on the 
sampling time and is subject to a large error, so only the 
"b' value is usually used. Values of *b* have been reported 
from 0.5 to 1.5, but it usually lies between 0.7 and 1.0 for 
tectonic regions. The "b* value observed at many lakes is 
given in Table 1.2 derived from the data given in Gupta and 
Rastogi (1976). It can be seen from this that the *b* value 
for after-shock sequences at lakes is about 1. They are on 
the high end of the normal range and usually higher than the 
'b* walue for the region. The foreshock ‘b*' values are 


larger than those of the after-shocks. 


The lake seismicity aiso has a large N/M, value where M 
is the magnitude of the mainshock and Mo the magnitude of 
the largest aftershock. N/M, values range from about 0.6- 
0.9. Gupta and Rastogi conclude the normal situation is one 
where large 'b* corresponds to ae small N/M ratio and 
Smaller ‘'b* with a larger N/M value. This is in contrast 
with the situation at reservoirs where large 'b* values are 


calculated with large M/M , (Papazachos (1974)). 


Little reliance should be placed on seismic parameters 
as they are difficult to determine accurately. The 
comparison of 'b* values with historic ‘b' values is always 
meaningless. For example, it seems impossible to compute ‘b* 


values for lcng periods while the means of detection are 
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changing. The 0.47 ‘'b* value for the Indian Peninsula has 
been taken from data covering 300 years and I fail to_ see 
its significance in a comparison to a "b* value computed by 
a modern seismic array. The Peninsula data was published by 


Gubin (1969). 


It is my conclusion that the calculation of *b* and 
a/o values is only useful in observing possible trends. 
Caution should always be taken in comparing seismic ‘lumped 
observables from one area with another, and then expecting 


the averages to apply to a third. 


The correlation between the water level and seismic 
frequency has been published by many authors (see Gupta and 
Rastogi (1976)). Hagiwara and Ohtake (1972) give typical 
results at Kurobe and show the highest correlation (0.41) 


exists, in their case, with a lag of less than one month. 


The stresses produced by changes of water pressure in 
the rock affect the total stress. Two contradictory models 
can be examined depending on the rock's water content. If 
the rock is dry the water will not affect the stresses until 
after it has diffused into the rock. In this situation the 
tock deforms elastically, with the water effects occurring 
at a later time. If the rock is saturated, the pressure 
changes are transmitted through the connecting pores over 
the area affected by the load. This initial pressure change 


is rapid in the neighbourhood of the source of the change, 
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The pressure gradients established cause diffusion which 
relaxes the gradients. The relaxation takes place after the 
loading, with the rate governed by the rock permeability. 
The rock is assumed to contain water in the work that 
follows . Since dams usually confine pre-existing rivers it 
is a reasonable assumption that there already existed ground 


water and that it extended to significant depths. 


Reservoirs are usually guite large, and may be tens of 
kilometers wide by several hundred long. For this reason the 
hydrologic regimes should be examined in an attempt to fully 
describe the conditions in the rock before the reservoir is 
filled. Large quantities of water may be transported by 
these regimes and the water pressure may differ 
Substantially from the hydrostatic head. In some areas the 
water pressure will be larger, and at others the pressure 
smaller than hydrostatic. These large scale features are not 
well known since they are usually not reported in the open 
literature on lake seismicity. Obviously hydrology is as 
important as the tectonic stresses, but since no data is 
available at the reservoir scale it has been ignored in the 
examples. It is possible to generalise the mathematical 


theory presented later to include these effects. 


One would assume the changes in seismicity would occur 
with the level changes and the later diffusion processes. at 
later times. Small foreshocks change the drainage pattern 


and hence the diffusion times required for fluid transfer so 
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correlations between lake level and seismicity may become 
confused at later times. The best correlations are usually 
seen during the first filling and several examples of this 
are shown in figure 1.10. The positions of most events are 
poorly known, and this particularly applies to the depth 
determination. At Vajont large scale slumping occurred, and 
Since the depth is unknown several processes may be causing 


the detected activity. 


It is not obvious from these diagraws whether the 
seismicity is related to the increased water lave or the 
rate of change of level. Both are related to the water level 
and, in diffusion processes, the rate of change of pressure 


is an important factor. 


Gough and Gough (1970) have indicated the desirability 
of using energy radiated, as well as frequency, in all 
correlations since this gives the correct weighting to the 
the larger events. In the examination of induced seismicity 
the strain released is proportional to the square root of 
the energy (Richter (1558), p368). The energy radiated is 
closer to the cause-effect mechanisp than are the 


readjustments shown by the seismic frequency. 


It will be shown later that it is possible to relate 
the stress produced in the rock by the lake with the water 
depth, the rate of change of pressure head, and the time the 


pressure head is maintained. The correlations between the 
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Figure 1.10: 


(a) Water level, frequency of events and 
energy radiated at Kariba are plotted as a 
function of time. The energy radiated is 
strongly dependent on the largest event 
within each of the four week periods. 
(Modified from Gough and Gough (1970b) ) 


(b) Histogram of the number of events and 
lake depth at Vajont in Italy. In October 
1963 a large mass of rock slid into the 
lake causing great damage. (Fron 
Galanopoulos 1966) 


(c) Seismic activity and water depth at 
Kremasta for 15 day periods. A _ sudden 
increase in activity is noted when the 
rate of filling was increased. (Fron 
Galanopoulos (1966)) 


(d) Flow hydrographs, water level and 
seismic frequency are shown at Koyna. The 
frequency increases about three months 
after the water level increased. In 1967, 
the pattern is confused, and on December 
10 a magnitude 6.5 event was detected. 
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seismic frequency and lake level shown in figure 1.10 are 


not always as obvious. 


The study of the seismicity is primarily concerned with 
changes in seismicity and Simpson (1976) has suggested the 
use of six categories. These are listed in table 1.3. It 
must be remembered that there are many times more lakes with 


no noticeable seismicity than are listed in table 1.3. 


Koyna, Kremasta, Hsinfengkiang, Kariba, Hoover and 
Marathon are placed in the category where there is a major 
change in the seismicity around the reservoir with a fully 
developed sequence of fore and after shocks. The case is 
quite clear at Koyna and Kariba where the N.O.A.A. Summary 
Tape and other research indicate that there were no events 


in the neighbourhood of the lake prior to filling. 


At Kremasta good correlations are observed between the 
smaller events and the water depth. However, as shown in 
figure 1.3, the epicentral magnitudes and distribution of 
the larger events before and after filling are not 
significantly different. The area was obviously close to 
failure and it is not clear whether the event was 
statistically likely to have occurred regardless of the 
presence of the lake. Oroville could be included in this 
category but the case for lake induced seismicity in this 


instance is less clear. Unlike the other cases where events 
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were observed soon after the initial filling, little 
activity was detected at Oroville until a few months before 
the major event. The time lag betwe2an the initial filling 
and the main event is the longest of the cases (about 7 
years). There may be scme connection between the downdraw 
just before the main event and its onset. The seismicity at 


Oroville is summarised in a paper by Bufe et. al. (1976}. 


(b) Minor Induced Earthguakes 


Simpson lists 12 cases where earthquakes of approximate 
magnitudes 3 to 5 have occurred and these are listed in 
table 1.3. In areas where prefilling activity was known the 
inclusion of Nurek, Monteynard, Bajura-Basta, Vajont and 
Benmore in the table 1.3 is based on increased activity 
coincident with the filling. At Talbingo and Kurobe 
prefilling instrumentation has shown that for 10 years prior 
to construction no earthquakes of the size measured after 


filling were detected. 


At the six reservoirs shown in table 1.3 instrument 
with sufficient sensitivity were available to detect an 
increase in activity below a magnitude of 2.5. A lack of 
instrumentation at other sites would indicate that this- may 


be a more widespread phencmenon than the data would suggest. 


(d) @ransient Changes in Seismicity 
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Simpson places five reservoirs in this category which 
have been associated with activity which began after 
impounding but which has since stopped or decreased. These 
are listed in table 1.3. Rothe (1970) has suggested biat 
hydraulic or crustal adjustments account for the seismicity 
at Oved Fodda and Camarillas. At Contra the activity stopped 
after the water level was lowered. The seismicity occurred 
after a rapid emptying and filling of the reservoir at 


Vouglans. 


These cases are very difficult to predict, but they 
indicate that the stresses were at a critical level and the 


region was close to failure before the lake was filled. 


As reported earlier, Jacob et. al. (1976) have 
indicated a decrease in activity at Tarbela reservoir during 
the initial stages of filling. This reservoir has only just 
been impounded so the change in seismicity may be a 
transient phencmemon related to the initial filling. This 
decrease is expected by the simple elastic modeis shown in 


figure 1.11. 


(£) Possible Cases of Reservoir Related Seismicity 


There are a number cf cases of increased activity near 
reservoirs where there is insufficient data to show a clear 


cause-effect relation. Simpson lists several in this 
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f 
category, however they are of little use in assessing the 


risk at a particular dam site. 


4 


The introduction of large amounts of water into the 
rock below the reservcir affects the stress pattern in a 
manner opposite to that due to the stress changes due to the 
weight alone. Laboratory studies of fracture and deformation 
of fluid filled reck (Handin et. al. (1963), Knutson and 
Bohor (1963), Nur and Byerlee (1971), Garg agi Nur (1973), 
Robinson (1959), and many others) have shown that the 
concept of effective stress proposed by ferzaghi (1923) and 
later supported by Hubbert and Rubey (1959), describes the 
behaviour of fluid filled rock. The effective stress is 
found by adding the fluid pressure, multiplied by a 
constant, to the applied stress. It must be remembered that 
pressure has the sign convention of an applied tensile load 


Since it tends to expand the solid. 


Oo = 


eff ~ °appliea t % P 


There is considerable discussion on the value of (Nur and 
Byerlee (1971), Garg and Nur (1973)) but for most rocks and 
engineering applications, acceptable results are obtained 
when the constant equals one. may be less for extremely 
high or low porosity and may depend on the stress. A 


description cf the variation of is given in section 4.1. 


The Mohr circle representation best describes how the 
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pressure affects the stability condition in different 
faulting environments. In these diagrams the normal stress 
is plotted against the shear stress with a circle drawn 
centred on the average of the maximum and minimum principal 
stress and radius equal to half the difference. The radius, 
relative to the principal stress axes, is called the maximum 
shear stress. Diagram 1.11a shows a Mohr circle with F the 


1 
largest principal stress and F., the smallest. The envelope 


3 
of Mohr circles at failure under moderate stresses has the 
form of a straight line 
T=C# on TAN © 
where tT , is the shear strength of the rock 


d? » is the angle of shear resistance 


o. , is the normal effective stress on the plane 
of failure 


C , is the apparent cohesion, and is the shear 
strength of the material under zero normal 
pressures 

C varies considerably and may vary from zero in a fractured 
sample to several hundred bars in an intact sample. In a 
large scale geological situation fractures would usually be 
present and C is small, and thus the rock has little or no 
tensile strength. Tan ¢ is the coefficient of friction along 
the failure plane. As angle » lies between 25° and 45° then 


the coefficient of friction is between 0.47 and 1.0 but is 


usually around 0.6 ( $=309). 


The failure envelope is found by triaxial testing of 


dry samples. Under extremely high confining pressures the 
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Figure 1.11: 


This illustrates the Mohr circle and 
Coulomb failure concepts. In the jacketed 
sample under compressive loads F, and F 
the Mohr circle (in the absence of fluid) 
is shcwn by the circle labelled ‘1'. 
Increasing the fluid pressure moves the 
failure circle to position '2'. When the 
Mohr circle touches the Coulomb failure 
envelope, failure will cccur along a plane 
at an angle 6 to F.. The effective 
stresses 07 103, are Shown by the final 
; E 
circle position, oes 
The largest principal stress is oO and 
the smallest o in each of tHe three 
fault conditions shown (b), (c), (d). The 
Mohr circles due to these principal 
stresses are labelled 1". Application of 
a vertical load, P, moves the circles to 
position "2" in a dry condition. The 
horizontal stresses are affected by 0.3P 
due to Poisson effects. If fluid is 
introduced into the faults at pressure P 
the circles move to the left to positions 
"3". The movement of the Mohr circles from 
initial to final positions are shown. In 
normal dip-slip and strike-slip faulting 
the area is made less stable by the 
application of the load and water 
pressure. Thrust fault regions are 
stabilised. 
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Shear stress exhibits a decreasing dependence upon the 
confining pressure as the failure mechanism changes fron 


intercrystalline slip to intracrystalline slip. 


Knowing the form of the envelope, the plane of failure 
can be predicted from the point the Mohr circle is tangent 
to it. Failure will occur on a plane © from the plane of F, 
in the two dimensional case. The angle 29 is measured from 
the axis of increasing © anticlockwise to the point the 


Mohr circle is tangent to the failure line. 


The introduction of fluid, water in the case of 
reservoir loading, into a dry sample under compressional 
loads PieF, , moves the Mohr circle by the amount of fluid 
pressure applied to the jacketed rock towards the failure 


envelope. 


The effect of a lake load on an area and the 
introduction of the water pressure head into different 
faulting envircnments is shown in figure 1.11. In each of 


the faults o is the largest and 03 the smallest of the 


ihe 
tectonic stresses. The Mchr circle before the lake is loaded 
is shown by circles labelled '1'. In normal faulting 9°, is 
vertical but is horizontal for strike slip and thrust 
faulting. The addition of the weight of the water at the 
surface increases the vertical stress by P=ogh. An increased 


vertical stress alsc increases the horizontal stresses by 


Poisson effects by an amount P\/(1-\) where V is Poisson"s 
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ratio. If vis a typical 0.25, then increasing the vertical 
stress by P increases the horizontal stresses by 0.3 P. The 
Circles '2* represent the combined effects of the tectonic 
stresses and the weight cf th2 water. If the water is 
introduced into the fauit region at the full pressure head, 
the stress circles are moved towards the failure envelope. 
The stresses are now effective stresses as the fluid affects 
all stresses equally. The final positions are shown by 


circles labelled *3'. 


The introduction of a fluid load onto a normal dip-slip 
fault or a strike slip fault moves the region closer to 
failure. Thrust faulting areas are stabilised by the same 
argument as the Mohr circle ‘'3" is smaller and thus is 


further from the failure envelope. 


The arguments presented here are for idealised 
Situations but the ccnclusions should apply to real 
faulting. In a realistic situation the effect of the load 
diminishes with the distance from the lake and the surface 
water head is not necessarily the same as that near the 
fault. These effects will be examined later. The arguments 
appear consistent with the observed faulting. Near 
reservoirs faulting has always been normal dip-slip or 


strike-slip except at Nurek and perhaps Tarbela. 


1.4 The Denver Earthquakes 


eaves ivatiet oud wee 
“agoeThe “fawea att a8 isdiae 
2 Pee e776 retro a 12) 


athenged igesod w O¢ao deol OF 
at SaeOio corps sat Key DN +e 
Sndic athe yd inet itiogfa oth! 


ai eedtd: ots «ei ider ai 9G 


Peoeilusds so8 4x 0 TE 
ieee or Yigan hlicde| ites ie 

heel aur 46 taetde aac abies Slapiioas é at 
reat ate Bos. gist ad3 sos) 
eae aaen weds ts ocak dee 
adieenpis alt asi van tan ee 
D0 sends tuna Devreade | 

se attest nin na>d | 


o 


ipneel ote i. 


8 


H 


- 
7 i 


; A 
4 2 
7 _ _ a e ' > a 


52 


The Denver disposal well of the Rocky Mountain Arsenal 
was the first clear example of earthquakes caused by the 
action of man. D.M. Evans (1966) called attention to a 
correlation between the vclume of fluid injected into the 
well and the number of earthquakes detected at a nearby 
Observatory. The events were up to magnitude 4 on the 
Richter scale and their proximity to Denver City was cause 
for concern. The well had been drilled through 3700 m of 
sedimentary rocks and finished in fractured Precambrian 
granitic gneiss. During the drilling it was found that the 
water pressure was at least sixty bars below normal 


hydrostatic head. 


In March 8, 1962, 1.6x10¢ m3 of fluid waste was 
injected into the well with a maximum injection pressure of 
about 38 bars. Earthquakes were detected shortly after 
pumping commenced. The well was then shut down from 
September 1963 until March 1964. Until March 1965 injection 
was maintained at essentially a hydrostatic head and a 
monthly average of 9.1x103 m3 was added. To increase the 
discharge, the pumping pressures were raised to 70 bars and 
pumping was suspended in 1966 owing to public concern. The 
earthquakes continued at a varying rate of from 4 to 71 
events per month with a magnitude 5.0, the largest of the 
series, occurring on 10 April 1967. From a plot of magnitude 
against the logarithm of the frequency this event was not 


unexpected since the relationship was stiil linear. On 9 
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August ancther large magnitude 5.3 event occurred, followed 
by a magnitude 5.1 event on 26 November. These events were 
Significant in that they were not expected from the linear 
relation between l1cg-frequency and magnitude obtained from 
the earlier recordings. This indicated a dramatic change in 
the pattern of activity. Since 1967 the number of events has 
decreased and no more large events have occurred. The well 
pressure, frequency and volume of waste injected are piotted 


as a function of time in figure 1.12. 


The U.S.G.S., the US Army Corps of Engineers and _ the 
Colorado School of Mines examined the seismicity near the 
Rocky Mountain Arsenal. A search of the historic records 
revealed no evidence of seismic activity before the weil 
injecticn in March 1962 similar to that which has occurred 
Since 1962. From a network of arrays of seismometers they 
concluded the epicentres lay within an elliptic zone 10 km 
by 3 km which encloses the disposal well. The focal 
mechanisms are consistent with right-lateral strike slip on 
vertical faults at N60W parallel to the major axis of the 


ellipse. 


Healey et. al. (1968) deduced from the mechanisms that 
the energy released must have been stored in the basement as 
elastic strain of tectonic origin. There are two reasons for 
this conclusion. Firstly, the magnitude-frequency relation 
is similar to that of southern California, which is a 


similar tectonic region, and secondly, because the radiation 
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Figure 1.12: The number of earthquakes per tonth and 
the well bottom pressure were plotted for 
Denver. The critical pressure determined 
from hydrofracturing tests is also shown. 
The seismic frequency increases when the 
pressure is above P critical. Between 1963 
and 1965 the number of earthguakes 
decreased significantly while the pressure 
was below the critical pressure. 
(Modified from Healy et al (1968)) 
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patterns suggest a zcne of vertical fractures existed with a 
northwest trend prior tc the injection of the fluid. These 


have now been observed in cored samples. 


A model was proposed by Healey et. al. (1968) to 
explain the faulting in terms of effective stress. The 
triggering of the region must have been a result of the 
fluid injection. In strike-slip faulting the largest and 
smallest principal stresses are horizontal. The pressure 
required to start fluid injection was found to be 362 bars. 
Above this value the injection rate increased remarkably so 
at 362 bars the cracks oren to fluid flow. One can conclude 
from this that the minimus principal stress, 0, , must have 
been 362 bars which is 93 bars above the initial fluid 
pressure in the well. The pore pressure before injection was 
calculated to be 269 bars by observing the rate of pressure 
fall-off with depth after pumping was ceased. The 
intermediate stress, On, 6 is due to the weight of fluid and 
rock. At 3671 metres, oe is 830 bars for saturated 
sedimentary reck of weight 0.226 batr/metre. The largest 
principal stress, O° is necessarily at least 830 bars. At 
the time faulting ‘was initiated the water pressure was 


measured at 389 bars. 


A ie O71 = 830 bars, og us 362 bars and P = 269 bars, 


the effective shear and ncrmal stresses on a plane at 60° to 


03 will be 7 = 203 bars and al 210 bars respectively. 


The angle of 60° was chosen since this is typically measured 
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in laboratory fracture tests. For a rock of frictional 
coefficient of 0.6, ee the cohesive strength would have to 


be at least 82 bars. 


During injection when the fluid pressure was 389 bars, 
x would again be 203 bars (independent of the fluid 
pressure) and ce 90 bars. The frictional tern on tan » is 
smaller by 69 bars indicating faulting would occur in a rock 
with a ge vaiue of 151 bars or less. Brace (1974) has shown 
the strength of rock decreases as the sample dimension 
increases. The fracture strength of quartz diorite decreases 
by a factor of five over two orders of magnitude change in 
the sample dimension while ccal's strength decreased by a 
factor of ten. The cohesive strength of sound rocks is 
around 500 bars so 150 bars may not be an unrealistic 


cohesive strength. 


An alternative to this model was also suggested by 
Healey et. al.. If the fractures already existed, then 
rather than failure occurring by the opening of new 
fractures, old cracks may have been reopened. Brace (1974) 
also discusses this question and points out that stresses 
and stress dreps associated with brittle fracture are 
generally much higher than those required for stick-slip 
faulting. Pre-existing faults nearly always occur in seismic 
areas and at Denver suitable fractures were observed in 


cored samples. 
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If 362 bars reopened the fractures, then the normal 
stress 0, on them was 362 bars. The least principal stress 
can be computed from the Mohr circle to be 

2 oe 0, (cos (2 +1) 
oa y= —_—— = 206 bars 
(1-cos (2q)) 
04 was 830 bars and © was again chosen to be 60°. If the 
fluid pressure was raised to 362 bars the Mohr circle is 
tensile for angles up to about 60°. Rocks have a low tensiie 
strength of only about 50 bars so pre-existing faults could 
easily open. Faulting would cease when the pressure has 
dropped sufficiently by permeating along the faults. This 
explanation of the failure was favoured by Healey et. al.. 
The results of the injection at Denver indicates a 
connection between fluid injection and induced seismicity. 
The later chapters of this thesis will try to examine the 
effect more quantitatively. The concept of effective stress 


has aiso been shown to be consistent with the failure 


observed. 


1.5 Phe Rangely Experiments 


The proximity of the Denver well tc populated areas 
made it impossible to conduct experiments on earth stresses 
there. An opportunity to do the experiments came at Rangely, 
Colorado where waterflooding of an oilfield has apparently 
triggered small earthquakes. A paper by Raleigh et. al. 


(1975) summarizes the results of the tests on fluid pressure 
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and the application of effective stress. 


The oilbearing Weber Sandstones are 350 m thick and lie 
at a depth of 1700 m. The Paleozoic sandstones and 
limestones lay beneath the Weber formation with the 
crystalline basement at 3000 m. The drainage pattern is 
aligned along a structure trending N30E. The Weber sandstone 
haS an average porosity of 12% and a permeability of 1 


millidarcy. 


Hydrofracturing tests were conducted in oil wells to 
determine the in-situ stresses. This technique is well 
described by Hubbert and Willis (1957) and Morgenstern 
(1962). The water pressure is increased until the tensile 
strength of the rock is exceeded in a packed off section of 
the well. At this critical pressure, the breakdown pressure, 
fluid flows into the rock and the pressure drops. The pump 
is then stopped and the fluid pressure drops until the pores 
close and the fluid pressure drop flattens. This is the 
instantaneous shut-in pressure and is equal to the normal 


stress acting across the fracture. 


At Rangely the least principal stress, Oz 6 equal to 
the shut-in pressure, was found to be 314 bars. fhe 


breakdown pressure, P was 328 bars and the initial water 


£" 
pressure, Po was measured to be 162 bars. The maximum 


principal stress, o, « can be found from:- 
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where fT is the tensile strength measured in the laboratory 
to be around 100 bars. From this formula o a= 552 bars. The 
intermediate vertical stress in this strike-slip region 
equals the overburden pressure of 427 bars. The coefficient 
of static friction for Weber sandstone is 0.81. The tectonic 
stresses were resolved along the plane of faulting 
determined by the focal solutions of nearby earthquakes to 
give shear and normal stresses of 72 bars and 342 bars 
respectively. The critical pressure required to initiate 


fracturing using these values is Po 257 bars... 


To test if the criticai pressure is correct the 
pressure at the bottom of four wells was raised to between 
235 and 275 bars between Cctcber 1969 and November 1970. 
There were 900 events in this period, 367 of them within 1 
kn of the wells. The wells were then sealed until 1971 and 
the pressure dropped in one hole from 275 to 203 bars. The 
number of events/month dropped from an average of 28 in the 
previous year to 1 event/month. In May 1971 injection was 
restarted. In well FEE-69, Ehegeeeesue: was raised to 265 
bars and the activity remained low at 1 event/month. Due to 
production charges at other parts of the oilfield the 
drainage pattern was altered at this time. A booster pump 
Was required to raise the pressure to about 275 bars. 
Between October 1972 and January 1973 with the pressure 
between the critical 257 tars and 275 bars, the monthly 


average rose to 6 events/month. Between January and March 
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1973 with the fressure abcut 280 bars, monthly average rose 
Sharply to 26 events/month. Backflow was initiated in May 
1973 and no events have occurred within one kilometer since 


that day. A summary cf the result is shown in figure 1.13. 


The consequence of the work described in this chapter 
are two-fold. It is clear that the presence of water 
certainly has an affect on induced seismicity and the 
effective stress concept has been shown to be sufficient to 
describe the triggering situation. It is on this basis that 
we approach the problem of seismicity at lakes. From the 
conclusions made it is necessary to propose a model which 
will explain the observed seismic features. The first model 
which will be introduced is one in which the weight of the 
water in the lake acts as a vertical surface load on an 
elastic layered half space. The second, and more important, 
is one in which the weight of the water acts on a permeable 
layered medium. The relation between these models and the 


available data is discussed. 
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Figure 1.13: The number of earthquakes per month and 


reservoir pressure are plotted at Rangely. 
The critical fluid pressure is shown. 
Notice a decrease in seismic activity when 
the well pressure was below the critical 
level. (Modified from Raleigh et. al. 
(1975) ) 
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The calculation of the response of a non-gravitating 
elastic half space to surface pressure is called the 
"Boussinesg" problem. It avoids the complications introduced 
by including the action of water in the rock and is a good 
starting point for the stress calculations which follow. The 
water is assumed to be a rassive surface load. It is treated 
as if it were an equivalent block of ice applying a vertical 
load to the surface of a perfectly elastic half space. The 
elastic behaviour is justified since, as will be derived 
later, deflections are cf the order of 10 cm over distances 
of 10 kms a strain of the order of 10-5 cm/fcom. This is well 


within the elastic limit of most rocks. 


Solutions of BoussSinesg froblems are common, and most 
engineering texts give examples, ( Timoshenko and Goodier 
1970) and a good treatment of geophysical problems was given 
by Farrell (1972). The application to lake loading was first 
done by Gough and Gough (1970 a,b) at Kariba. They 
represented the lake as 1302 point forces equal to the 
average pressure over rectangles 2.2x2.3 km. A considerable 


effort was necessary tc determine these forces and a 
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description of attempts to find an easier method is given in 


a later section. 


The stresses and displacements for a point force 


-F(0,0,0) at a point P(x,y,z) are given by: 


he! mee 2 
ole = Og Sin 6 + ox cos 6 
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and the vertical deflection: 
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These are derived in Timoshenko and Goodier (1970) and are 
used by Gough and Gough. In these equations: 


R2 = x2ty2+z2 
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ce = x2+y2 
tanod = y/x 
Vv Poisson's Katio 
Y Young's Modulus 


Kariba reservoir is about 200 km long, 25 km wide and 
is 80 metres deep. The largest vertical deflection at 3 km 
below the lake was 23 cm. The largest vertical stress, 06;, 
was about 8 tars under the deepest section of the lake. fhe 
maximum shear stress was about 2.1 bars again under the 
deepest part of the lake. fhere is little point in 
reproducing contours of these results at this stage. Later 
an alternate method will be developed for finding stress and 


displacement . 


The elastic constants were taken from values published 
by Birch (1966). The values were v=0.27 and Y=0.85 Mbar 
which are consistent with a sound, intact rock. The 
displacements at Kariba were compared with road levelling 
done near the dam wall (Sleigh 1976). The results are shown 
in figure 2.1. They show an excellent agreement between the 
calculated and observed vertical deflection. The agreement 
confirms that analytic techniques can be applied to rock 
deflection over a scale of 10's of kilometers. The 
relevelling has been repeated by Sleigh (1976) and little 
change is seen. Some deviation is noted in areas of Karroo 


sediment ; this is probably due to hydration. 
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Figure 2.1: 


Observed and calculated deflection 
measured alcng a roadbed near Kariba in 
Rhodesia. The calculated deflection was 
computed by the Gough and Gough (1970a) 
point force technique. 

(Modified from Sleigh (1976) ) 
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The solutions obtained by summing the contribution of 
point forces has several advantages as well as 
disadvantages. As Gough and Gough showed in their 1976 
analysis of Cabora Bassa, the point forces need not be 
equispaced. In sections of the lake which are long and thin, 
or, in which the depth changes rapidly, higher sampling 
rates may be made. The solution for the displacement (2.1) 
is of a 1/R form and is numerically indeterminate close to 
the point force. Gough and Gough (1970a) avoid these 
Singularities by computing no closer than 4/3 of the force 
separation. This is also justified since the approximation 
of the lake by point fcrces breaks down closer than this. 
The assumpticn of a flat surface and vertical loading is 
violated too close to the surface. The largest disadvantage 
of the technique, acknowledged by Gough and Gough, is the 
cost of computation. 14 minutes was required for a vertical 


section on an IBM 360-67. 


The technique of summing point force contributions has 
been applied by Beck (1975) to Oroville in California. Fron 
the bathymetry, displacements and maximum shear stresses 
were conputed. The largest shear stress was 3.4 bars under 
the deepest part of the lake. This value is larger than 
Kariba's maximum shear stress of 2.1 bars since Oroville is 
about twice as deep as Kariba. The maximum vertical 
displacement was estimated to be around 5.5 cm, again at the 


deepest lake section. Beck has extended the point force 
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soluticn to the surface by apparently trying to compute 
between the forces. He admits the deflection contours are 
seriously distorted by local effects. These distortions have 
been smoothed from the final results. These distortions 
could have been predicted during computation due to the 1/R 


Singularity in displacement. 


Another technique has been developed by Lee (1972). The 
Boussinesq problem was solved by a Fourier-Bessel expansion 
in cylindrical coaqrdinates. The lake was approximated by a 
pattern of elemental loads, each shaped as sections of an 
annulus. The method is shown to be able to compute the 
displacements with the elastic parameters varying with 
depth. The expression for the vertical displacement is: 

Be 


4 


5 T, (1+2v) (1-v) a : 
Tay) 1 (V(b.)-V(a.)) + 
elements,i ies = + 


A. 
ak 
+ 2°(W(b,)-W(a,)) Jao 
i on 
where: 
a 
V(x) = [x? + 7 - 2rx cos(6-a) + 2°]% + ot COs Gad) x 
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x ln Co ae x (2cte +x? 20x cos (8-a) ) 7 
ie ray 5 ATE yee Se ae J 
(Zero Sa ne (b=) ) 2dr Isa npe(Oso) 
and; 
ian [xr cos(6-a) - proc’ 
Ee EI PRED Te) ale 
fiz -+r% isin~(6-a))[2etre=2rx Cos (0-d) 2%) 


V(t, 0,z): is the vertical displacement at (Lr, 9,2) 
relative to an established cylindrical coordinate 


origin. 
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T; is the average pressure over an area bounded by 
radii a andb and subtended between the angles A and 


B . 


Y: is Young's Modulus and y is Poisson's ratio. 
The summation is taken over each of the elements used in the 


approximation of the lake load. 


Lee did the computation by a five point Newton-Coates 
technique for Lake Gordon in Tasmania. The lake's shape was 
approximated by 51 elements and the computation was done at 
598 points. No description is given of the treatment of 


Singularities or the efficiency of the computation. 


2.2 A Transform Solution of Boussinesg's Problen 


For any material the conservation of forces and noments 
leads to the equilibrium equation: 
ees ei. 5, 
- 13,7) 
where o,, is the stress ccmponent applied in the direction i 
13 


= a0 ie eee 


to the plane whose normal is in the j direction. F. are the 
body forces in the i direction. Normal index notation will 
be used with indices having values 1,2 or 3 unless otherwise 
specified. Repeated indices indicate summation over the 
index, and commas indicate differentiation with respect to 
the variable. There are six independent stresses due to the 


symmetry of the stress tensor, SO: 
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For small deformations the material behaves elastically with 


stress linearly related tc strain. 


O.. = ie 
‘ia ea, iki ae 
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oa 2 2 ee 4 Ox ) 
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where u. is the deformation in the i direction relative to a 


cartesian cocrdinate system u(x ,, X., X,)- 


The a are constants of proportionality. The 21 
independent constants can ove reduced to 9 if the material 
has symmetry of properties around 3 mutually perpendicular 
axes. If the elastic froperties are the same in all 
directions (isotropic) the number of elastic constants 
reduces to 2. In all the following work Poisson's ratio, v, 
and Young's Modulus, Y, are taken as the measured elastic 


properties. For an isotropic elastic material the stress- 


Strain reiaticnship reduces to: 


x Vv 
rf. a ee ee 
dl | Cy) Fti=20) rie Spyegs rae re 


To examine the stresses created by the application of a 


surface load to an elastic half space the equilibriun 


condition is differentiated to yield: 

V(Veu) + (1-2v)V2u = 0 ees 
Landau and Lifshitz (1957 pp26-29) show that this equation 
can be reduced ty replacing the displacement vector, u, as: 
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expressed in terms of another harmonic vector: 


0g Pp) 
x OZ y dz 
2 
= Z 
u Fy 0 V Jy = 0 


Substituting into (2.2) yields: 


ag dg 
ee mae, xX 
= f . my = 
ci-wy 2 aoe oy ee eB 


where ww is ancther harmonic scalar. A Cartesian coordinate 
system has been chosen with z positive below the surface on 
which the forces are applied. For the ccmplete description 


and ~ are required. These are found 


of the solution Tx@Iyet, 


from the bcundary conditicns at z=0. At the surface Oe et 2 


where P; are the applied forces in the x, y, andz 
directions. The negative sign is the normal convention to 
indicate that the surface is under compression. Tensional 


forces are positive. 


Expressing u in terms of the four unknowns, and 


substituting into (2.2) yields the conditions at z=0: 
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ag 2(1+v)P 
cs wpiecies oY 5 z . 
Gat ae Se -2.6 


Z=0 
A further condition is required to uniquely define the four 


unknowns. Landau and Lifshitz show that it is convenient and 


noncontradictory to set: 


0g 0g 
5 ty ee a ak eae, 
(ek: 2v)£, = a y ) + 4(1-v) ao eo =2 27 


The first two relaticns of (2.6) become: 


2 
3 J, — 2(1+v)P, 
2 > Y 
See z=0 
-2.8 
2 : 
ag 2(1+v) P 
ee é ce ey eo ok ee a 
Z Y 
OZ 
z=0 


Green's theorem is used for the solution of y and tf. The 
theorem may be simply stated. If a function is harmonic, is 
zero at infinity, and has derivative df/dz at z=0, then: 


AGrtiae) = ri tea y! Ze) a 


z=0 

The integral is taken over the area of definition of f. R’is 
defined as ( (x-x')2 + (y-y*)2 + z2)  . Treating the load as 
a vertical force, Poe or a free surface with Set oe and 
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Lil lets ow 
From the a rae Of 37 « jx and jy can be obtained by 
integrating st over the variable z and then 


differentiating. 


The solutions used by Gough and Gough (1970a) for the 
point force are easily obtainable at this stage by examining 
the z displacement due to a point force F.6(x").6(y") From 
equation (2.4): 


Up eer & 4 4+ wD 
Z Z Z 


ee =4) fa Zz an po 
4(l-v) “z 4(l-v) 902 OZ 
ie and av can be found from (2.9). 
OW 2). (itv) (ee) ) ee 
OZ 2m¥ ie 


22 Clee ey 


vA oil TY r 
ee 2il=v he 
OZ TY 
Therefore: 2 2 
ee ee 
vA 2nY ie 3 


r 


This agrees with equation (2.1). Similar techniques are used 
in the determination of u ,v . The integral (2.271-4, p86) 


of Gradshteyn and Ryzhik (1965) is necessary for these 


computations. 
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The 45 stresses are found by differentiating. 
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For the lake , the surface load is continuous. Equation 
(2.9) can be solved using a Fourier Transform to remove the 
integration variables x*',y'. Multiplying both sides of (2.9) 
of f. by the transform variable and integrating over x,y 


gives: 
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p and q are the Fourier transform variables, or wave 
numbers, for x and y. Using the substitutions: 
et) B= eCos y eee ees LG 


p = u cos® gq = u sind 


and the notation: 


io 
Cpe ay 2) = ff aeryrarel? e Gy dxdy 
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the transformed equations become: 
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Using Gradshteyn and Ryzhik (henceforth G.R.) (3.915, p482), 


(6.554, p682) 
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sWsxerey= ; 
ax = ie T J&B) (J, is a Bessel) 


fONCLION 
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Solutions for f. and W can ke obtained: 
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By using the transform of (2.4) and (2.5), the equation for 
the displacements can be found. 
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Where k2=p2+q2. To compute the stresses, derivatives with 
respect to x,y and z are required. The z derivatives can all 
be done analytically, however the x and y derivatives are 


‘simplified by a property of Fourier transforms. 
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The stresses and displacements can now be computed at a 
depth z. If we know the load P(x,y) the P(p,q) can be obtain 
by the fast Fourier transform (F.F.T) algorithm presented by 
Cooley and Tukey (1965). This aigorithm is available at most 


large computing centres. 


A singularity appears in the determination of the 
displacements. This difficulty exists for p=q=0. P(0,0) 
corresponds to the integral of the load over the area 
transformed. It is the average load and in finite transforms 
is always non-zero. In order to avoid such computational 
Singularities, we either subtract the D.C. level from the 
entire grid P(x,y), before transforming, or do the inverse 
transforms in cylindrical coordinates. The first device is 
forced on us by cur desire to use a F.F.T. algorithm. If the 
inverse transform is carried out in cylindrical co- 
ordinates, Par. ts methods are not applicable. After 


inversion the D.C. contribution is added to the final 
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solution for the displacement. 


The deflection and stresses due to a uniform square 
load can be derived from results in Love (1929). Love's 
solution can be computed at the surface if the edges of the 
square are avcided. The total solution can be found if the 
bounding area of the square is large relative to the size of 
the lake load. This shculd always be the case since it 
reduces to magnitude of the load on the square. For a 
rectangle of sides 2a ky 2b centred at the origin of (x,y) 
each of the functions are evaluated four times due to the 


contribution of each edge of the rectangle. 


f(x,y) =f (xta,y+b)-£ (x+a,y-b)+f (x-a, y—-b) -f (x-a, yt+b) 
The displacements can be shown to be: 
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The stresses are easily derivable from these: 
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ee a 
Oy = aa ((1-2v) 


The potentials y and V are: 


x = { frcstz+n ax dys 
¥ -[f dx'dy' 
R 


where R2= (x-x!') 2+ (y-y!*) 2+z2 
V may be shown from G.R. (2.261 p81, 2.266 p84, 2.267-1 p84) 


and partial integration to be 
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The remaining derivatives of V, can be calculated in a 


straightforward manner. 


Figure 2.2 shows the vertical displacements below 1 bar 
pressure on a square of 284.16 km. This length had been 
chosen to agree with the Kariba results presented later. 
Figure 2.2 shows the ccmputation at O km and at 13 km deep 
are only slightly different. This indicates the displacement 
changes slightly close to the surface. Figure 2.3 shows the 
stress 10 Km below a rectangle on which is applied a one bar 
pressure. Figures 2.4 and 2.5 show displacements and 
stresses below a sguare with a side of 100 Km. The same 
elastic parameters used by Gough and Gough (1970a) were 
chosen. The bar figure in the lower right of each of these 
diagrams is one tenth the lengh of the square. the length of 


the bar is given in kilometers. 
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Figure 2.2: 


Vertical displacements at several depths 
below a 1 bar pressure applied over a 
Square of side length 284 km. The elastic 
parameters of the halfspace were Poisson's 
ratio of 0.27 and Young's Modulus 0.85 
Mbar. The equations derived by Love were 
used. 
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Figure 2.3: Stresses at 13 km below a unit pressure 
applied to a square of side 284 ka. 


Figure 2.4: Vertical displacement at several depths 
below a unit pressure of 1 bar applied to 
a Square with sides of 100 kn. 
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Figure 2.5: 


Stresses at 10 km below a square of side 
100 km on which a unit pressure of 1 bar 
is applied. The solutions were obtained by 
the equations derived by Love. 


She2 2e0BARS.) 


DEPTH= 10 WITH SIGMA=0.27 YOUNG=0-.85MBAR 


SPGAY CBA Oo) 


DEPTH= 10 WITH SIGMA=0.27 YOUNG=0-85MBAR 


SIGXX (BARS) 


DEPTH= 10 WITH SIGMA=0.27 YOUNG=0-8S5MBAR 


ee ey 


SEND nal aes ina 


DEPTH= 10 WITH SIGMA=0.27 YOUNG-0.85MBAR 


9IGXZ (BARS) 


DEPTH= 10 WITH SIGMA=0.27 YOUNG=0.8SMBAR 


OI PYZ, 1 BARS 


DEPTH= 10 WITH SIGMA=0.27 YOUNG=0.8SMBAR 


tp: Tk 
eae. GOT eats in ic 


oe ——— <i apegey 
—- : 
ele ed oh aa a 


—* 


Naan eon 


89 


2-3 An example of Transform Techniques 


The computation of the stress and displacement was 
undertaken for Lake Kariba. The point force magnitudes and 
positions obtained from the bathymetry were supplied to us 
by Dr. and Mrs. Gough. The pressures were the average depths 
over rectangles of 2.2x2.3 km. Over 1000 point forces were 
used. This data was placed onto a regular grid of up to 
256x256 nodes. The 'D.c." level of the load was removed and 
a two dimensional transform was taken using a library F.F.T. 
routine. The F.F.f. programme was a standard IBM routine 
"HARM® supplied in their Scientific Subroutine Package. 
Singleton (1969) has published a listing and explanation of 
another F.F.T. programme. The transformed load was 
multiplied by the appropriate functions and the inverse 
transform was taken. A grid of 64x64 element was usuaily 
found adequate although larger grids can be handled 
efficiently. On the IBM 360/67 only ae single stress 
component or displacement could be computed at one time for 
a 256x256 grid. The difficulty involved the storage of the 


65536 source point forces. 


On the Kariba example we found the results using the 
F.F.T. routines cost about 1/8 of what the Gough method cost 
(64x64 grid). If the accuracy criteria are relaxed to 5 per 
cent agreement, this fraction can become as small as 1/20 
(32x32 grid). The loss of agreement is directly related to 


how frequently the lake is sampled by the grid nodes. There 
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are penalties when the grid is reduced. Transforming 
digitised data imprecisely on a limited grid may introduce 
aliased components of displacement and stress near sharp 
gradients of load, (see Kanasewich (1976)). It should be 
emphasised that routine processing and improved data can 


reduce these errors. 


The displacements and eigenvalues of the stress matrix 
were computed at several depths. From the eigenvalues’ the 


Maximum shear stress was calculated: 


where O and ee) are the maximum and minimunm 
max min 

eigenvalues of the stress matrix. The displacements, maximun 

shear stresses and shear planes at 0.,3., 13., and 30., Ka 

are shown in figures 2.6, 2.7, 2.8. The derived stresses are 

everywhere within several per cent of the results of Gough 


and Gough (1970a). Since more points were computed for the 


contours the extra detail is probably real. 


The deflection at 0. km, shown in ¢igure 2.6, differs 
little from that of Gough and Gough at 3 km. This result 
reinforces the comparison between the computed and measured 


relevelling data shown in figure 2.1. 


2-4 Generalisation to Regular Loads 


The prime advantage of the transform technique is its 
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Figure 2.6: 


The naxinuag shear stress, vertical 
deflection and the planes along which the 
maximum shear stress acts are shown at 13 
km below Lake Kariba. fhe position and 
bathymetry of the Lake are also shown. The 
two planes are each described by the 
azimuth of the symbol and tke dip angles. 
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figure 2.7: Vertical deflections at several depths 
below Lake Karika. The contours at 0 kao 
depth are at 2 cm intervals. 
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Figure 2.8: The maximum shear stress at several depths 
below Lake Kariba 
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ability to handle irregular shaped loads efficiently. In 
some circumstances a large computer may not be available and 
an analytic treatment is desired. The integrals involved are 
usually difficult tc solve and little attention has been 


given to the method. 


in some instances the reservoir may be approximated by 


one or more ellipsoids. The ellipsoid: 


x-x', 2 -v! 
(FA + CEL) 2 4 (2% 2 1 


would transform to give the fressure: 
2 
PG abC ee Tr 7372”) 


6 
where §2=(ap)2 + (bq)2, p is the load density and J 5/18 the 


Bessel function of order 3/2. A program was written to fit a 
small group of best fitting ellipsoids for Lake Kariba by 
least-squares methods. It was slow. For the most general 
ellipsoidal crientation the fitting alone required almost 
the same computing time as the complete solution by 
transform methods. This penalty may be acceptable in 
computers with small memories in which large transforms 


cannot be done. 


The technique is alsc amenable to two dimensions. If 
the lake is straight, with the length about three or more 
times greater than the width, two dimensional analysis gives 
excellent results. The equations are easily obtained by 


setting g=0 in (2.12), (2-13)- 


<1 ola ov euotes 
bisa ee aa on tort “a es 
odd ah, > Bee Yi tee Yaek sit 2k) ASG) * Stgnn==0 ox0du = 
6 +i2 of nSdteiaw e5u wsipesg a -Sye webs ie aotspar? feasoe . L 
gi 6S6digsA ede, 205 abpsorgi tie pol sgk% send Lo geozp Lowe vi 
fpxguey 208 Bed AE. .wele Bse “I -aboddon aetsupe-sesal * 
taumls be2enper SielLs BRERSES SAD ‘Gubtes as kx ae 
1 Hoidetom s2eiyKoo 4% 2s cube Bikswabo> vase ple feo 
gi Sadedusods od yon (tisiey aat bhodden rected 
eetctensss Spast doidv oi arinowew Bidwe dsie ier oH 


-saob od ¥ i 
Br veopceneais Govt of afdenvae 
enih Soimenie Secs atvodt sig dele Boars | 
Deiat is cint ae aw one 
Yd Benkssti> ybicse oe acoktoupe ait .attvces $aeLteone 
+ ters) “es a hi Ger 


gg 


98 


2-2 Generalisation te a Layered Half Space 


se ee ee ae ee 


The generalisation to the layered half-space follows 
directly from the derivations in section 2.2. The extension 
is made by considering the boundary conditions at each 
layer. Equations (2.2) to (2.8) are applicable, but as f.. 


and fy can nc lcnger ke neglected at each interface equation 


(2.9) no longer applies. 


Each of the four unknowns ), foe Je Jy are transformed 


and are separated into two exponential components as 


follows: 
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From equation 2.7 it is easy to obtain: 
(da) 


fe z (pg + qd=) 
es at Tok 2 ewe oe Ua ey 


where the v form is assumed to follow directly from the 
definitions above. The p andq are the Fourier transform 
variables in x and y with k2=p2+q?. The displacements can be 
obtained from the equation (2.4). The results can be 


summarized in matrix form: 
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The complete solution for the displacements in matrix forn 


is: 
is +A Seed coi ae 
Per LU) Bowe dea me unn atic ieeoree 
u 
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u 
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where ; 
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Dy cerwa pepe is 
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The expressions for the stresses can be obtained from the 
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equation (2.2). Using the matrices [L] above, and the 
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differentiaticn properties Fourier transformed variables, a 


Similar matrix notation may be obtained for the stresses. 


x wat = = 
Op, 5 = Gh e kz2 a 5. l akZ 
ak Jb 7 1j 


For convenience in the introduction of the boundary 


condition the six stress components are separated as 


follows: 
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i Y 
LSA = a Geen ATES) 


-ip[p* (zk+1)+2k?(v-2) ] ~iqlp? (zk+1) #2vk7] [p? (zk¥(1-2v) )¥ 
2vq"] 

2 iL a ze 

(v-2) ] [q” (zk+(1-2v) ) + 


2vp*] 


~ip[q? (zk+1) ¥2vk7] ~iqlq? (zk+1)+2k 
: 2 _ 2 ; 2 — 2 — 
-iq[p° (zk+1)+2(1-v)k*] -ip[g (zkhat) 2-9) k= [pal zk 25) 7 


The layered space identified in the figure 2.9 has N layers 
each with Young's Modulus and Poisson's Ratio Y ov and 
depth at the top Zoe The N'th layer is a haif space. fhe 
matrix solutions for displacements and stresses at any depth 


Z, in layer i, are written below. 
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computation is in the i layer at depth z below the free 


surface. The 1 is the 3x3 unitary matrix and 0 is the zero 


3x3 matrix. 


At each boundary there are 6 boundary conditions. The 3 
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Figure 2.9: 


Shows the configuration used for 
calculating the displacements and stresses 
in a layered half space. There are N-1 
layers over a half space each with Young's 
modulus Y and Foisson's Ratio - The top 
of the layer iis at z with z being the 
free surface on which the load is applied. 
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displacements are continuous if the layers are to stay 
welded together. The three stresses o , 060 , O must be 
ZZ XZ yz 
continuous across each interface to ensure no acceleration 
of the interface. The top and bottom interfaces have special 
conditions imposed on them. The surface is free, so we 
conclude o =o =0. The o is the applied ioad at the 
XZ yz Ze, 
surface. Below the lowest interface lies a half space. To 


ensure the variables are bounded with increasing z, the 


ig]. are necessarily zero. 


The boundary conditions allow the determination of 


(917-18); in all the layers: 


kez. =a, 
oe —, 
g* et: 0 M M 
7 =kz 
= oe _ 
g 0 e L L 
in—s]s i-1,z, 
-kKZ 
- a+ 
mM M e 0 g 
+ L “kz, t 
L L 0 e g 


The [s], =0 and foie can ke computed by chain application of 
win 
the identity above to the top surface. Knowing Ig]. the 


displacement and stresses can be calculated everywhere. 


The ccmplex matrices [L] , [M],- and [N] can be 
separated into a real matrix, pre and post multiplied by 
complex diagcnal matrices. The inversions can be found using 
real algebra routines. Caution must be taken to include the 


complex matrices when solving for the [o] matrix as this is 
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complex since the boundary conditions, P , are also complex. 


| Analytic inversions were attempted on the matrices 
L, M,N. These proved difficult to obtain and little 
Simplification was achieved. In ccmputation it was found 
more convenient to compute each of the matrices L,M,N and do 
the inversion numerically. The results obtained agreed well 
with the elastic case when the three layers were chosen to 
be the same. This is illustrated in figure 2.10, where the 
contribution of the square load has not been included. 
Figure 2.11 shows the results when three different layers 
were chosen. The technique was found to be unsatisfactory 


due to numeric instability in most of the models chosen. 


A computation difficulty arises again for the zero 
wavelength compcnent. Love's solution may be adapted into 
the layered mcdel to remove the singularity if the average 
pressure over the grid is smail. For a square much larger 
than the lake little difficulty will be experienced in the 
interpretaticn of the results. The solutions in the region 


of interest will only be slightly in error. 


In this chapter the Boussinesq problem has been 
examined for elastic halfspace. The Fourier technique 
developed has been shown to be more efficient 
computationally and is able to be calculated at the surface. 
The theory was extended for the elastic layered haif-space 


using Haskell-Thompson matrix techniques. 
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Figure 2.10: 


The vertical deflection is shown for the 
layered elastic halfspace (a). This agrees 
well with the solution fcund for the 
eguivalent halfspace without interfaces. 
The effect of the average pressure has not 
been inciuded in these diagrams so they 
differ in magnitude from the results shown 
in figure 2.7 
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Figure 2.11: Vertical deflection below Lake Kariba 
found using the 3 layered half space with 
the parameters shown. The 'DC' MJlevel has 
not been included in this case. 
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The model does not account for the effects of water, OL 
allow for the concept cf effective stress . These will be 


included by using the work of Biot (1940). 


The loading of reservoirs changes the pressure of the 
water contained in the rocks. Since reservoirs are built on 
rivers, groundwater is usually present and is observed by 
the presence of a4 watertable. Tne weight of the reservoir 
compresses the pore volumes and this effect must be included 
in any model of loading processes. The reservoir also allows 
water to enter the grcecund by leakage from the bottom raising 
the groundwater pressure. The modelling of a leaky reservoir 
on the porous elastic halfspace can be described by the 


concepts of effective stress and diffusion. 


3-1 Biot"s Consolidation Equations 


Biot (1941 a,b,c) developed the basic equations 
governing the behaviour of solids containing fluids. The 
theory was developed for the study of consolidation of 
foundations in clay and sandy material but is sufficiently 
general that it can be applied to rock. Compaction of soils 
is described in articles -y Hedberg (1936), Skempton (1970), 
and Marsal and Philipp (1$70). A recent review article (Rice 


and Cleary 1976 a,b), concludes that recent work has not 
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greatly improved on the classic (Biot 1940a) and later works 
of Biot (1955, 1956a, 1973). This is the analysis which will 
be followed here. The extensive engineering literature 
concentrates on. two-dimensional plane strain loads, and on 
problems with axial symmetry in three dimensions. Some of 
the recent applications inevader that by Nur and Booker 
(1972) and Booker (1974) to fault relaxation and by McNamee 
and Gibscn tc settlement (1960a, 1960b, 1957). A description 
of the solution cf engineering problems by finite elements 
will be given ina later section, but it is apparent that 
very few analytic solutions exist. Biot's theory will be 
applied here to irregular loads which change with time in 


both shape and magnitude. 


Hooke's equation of elasticity in the absence of fluid 


may be written as: 
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v is Poisson's Ratio and Y is Young's modulus. The water 
pressure does not alter the shear components of stress but 
affects the rest equally. If o is the water pressure, and H 


a parameter describing the effect cf the pressure. 
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The sign conventicn used requires the water pressure to be 


positive with increased pressure and the stresses negative 
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in compression. Three assumptions have been made -- 
(2) the material is isotropic 
(ii) a linear stress-strain relation exists 


(111) the strains are small. 


The amount of water in an element of volume is related 
to the six applied stresses and the fluid pressure. If the 
material is isctropic then the shear stresses cannot affect 
the water content. As there is symmetry of the material 
properties the contributicn of each of the remaining 
stresses must be equal. In the first order theory Biot uses 
a linear dependence of the water content, @ , on the 


pressure and stresses. 
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8 is defined as the increment of water volume per unit 
volume of the material. H, and R are two physical constants 
of the rock. If it is assumed that consolidation is a 
reversible process, it is possible tc show that H,=H by 
energy conservation arguments. The constant R allows for the 


inclusion of air with the water in the rock skeleton. 


Equation (3.2) and (3.3) can be written in a more 


useful way as: 
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The four material constants G, v, H and R and the derived 
constants a and Q need to be physically interpreted. G and 
vo are the equivalent elastic moduli of the solid skeleton in 
the absence of excess fressure. This condition can be 
established if experiments are conducted at a sufficiently 
Slow rate to allow excess fluid to diffuse out of the rock. 
In order to interpret H and § consider the sample to be 
enclosed in a rubber jacket. In the absence of external 


stresses and with the pressure in the jacket reduced by -o : 


where 6 is the amount of fluid removed and ¢ is the change 
in the volume of the solid. The constants Q and a are 


meaningful when the solid is not saturated. 


The constant a is the ratio of water volume squeezed 
out to the volume change of the solid if the latter is 
compressed while allcwing the water to escape. The a of a 
saturated material is equal to 1., and for dry material it 
is 0. Irobe et. al. (1974) use a values between 0.83 and 
0.98 in a description of the behavior cf soils. Jaeger and 


Cook (1976) indicate that aq has a value of one for rocks. 


1/0 is the measure of the amount of water which can be 


forced into the solid under pressure while the volume is 
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kept constant. A standard soil test affords a method of 
determining Q. A compressive load is applied through a 
porous slab to a laterally ccnfined column of soil so that 
the water is allowed to drain away. The ratio of the axial 
strain to the applied load is equal to the “final 
compressibility coefficient". This is shown by Biot (1941a) 
to be 
_ _(1-2v) 

2G (1-v) 
If the test is repeated with the strain being measured 
before the water has had time to diffuse out of the solid, 
then the ratio of the strain to the applied load equals the 


"instantaneous compressibility". 
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Q can be determined from these conditions with a knowledge 
of - Biot has shown that the initial compressive strain is 
large compared to the final compressive strain for soils. 
The ratio cf instantaneous to final compressive strains is 
typically 0.01. On this basis Q is large, and several orders 
of magnitude larger than Young's modulus. Jaeger and Cook 
(1976) suggest that as the compressibility of rock-forming 
Mineral is of the order of 0.1% per Kilobar, it is 
reasonable to neglect the volume change of minerals at 
laboratory pressures so that a Q value of infinity is 


suitable for rocks. 


Rice and Cleary (1976a) modified the choice of derived 
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constants so they more closely resembled bulk modulii and 
have introduced the effective stress law of Nur and Byerlee 
(1971) rather than that of Terzaghi. The lack of information 
an any of the constants cver scales of the size of lakes 
(10's of kilometers) is the largest single weakness of 
theoretical analysis. In the work of Rice and Cleary (1976b) 
on faults, the parameters may be monitored more closely, 
however, for jlake analysis the simplest macro-parameters 
will be used. The notation of Biot will therefore be used 


throughout. 


Equation (3.4) can be differentiated and simplified to 
yield the first of the basic equations required. 
ena us + .(2n-y)e iad - ach. = 0 -anG 
It 1s assumed the water's behaviour is governed by DArcy's 
Law where the fluid flow velocity is linearly proportional 
to the fluid pressure gradient. For an isotropic solid with 


coefficient of permeability K 


2 ~3.7 

The assumption must be made that the flow is not turbulent, 
but its wide use in oil and hydrological studies implies 
this assumption is justified. If the water is 
incompressible, the increase of water content in an element 


of the solid equais the flow into it: 
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Using (3.5), (3.7) and (3.8) 


~3.9 
The Biot Ccnsolidation Equations (3.6) and (3.9) fully 
describe the deformation of a fluid-filled elastic material 


due to an applied load. 
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The equations (3.6) and (3.9) can be solved using the 
displacement functicns of McNamee and Gibson (1960 a,b). 
Other solutions can te applied (see Rice and Cleary (1976a)) 
but are generally limited to two dimensions. Taking the 
divergence of the vector equation (3.6): 

V* (2nGe - ao) = 0 
If S is a harmonic functicn ( Vy2S=0) then, 
2Gne - ao = -2G = 
It can be shown that ancther scalar function E can be 


introduced which will satisfy this equation if 
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Substitution of (3.10) into (3.9) and time scaling by t*=tC 


leads to the two differential equations 
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The time scale ‘C* is linear in the permeability, K, so one 
result may be examined for several permeabilities. For 
example the results computed at 1 year for a rock of 20 
Millidarcy (20x10-!%cm/s) will be the same as a rock of 2 


millidarcy computed at 10 years. 


The equations (3.11), (3-12) are solved by Fourier 
transforming ia the xX and y variables and Laplace 


transforming in 't*. The notation used is: 
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With k2=p2+#g2 (3.11) and (3.12) reduce to: 
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The condition at t=C corresponds to no initial 
compressibility and no excess pressure. This assumption is 
made on the kasis that the results will reflect the effects 
of the ancmalous load alone. If the initial stresses, fluid 
pressures and flow velocities were available they couid be 
included at this stage. Body forces should be considered in 
this case as they will cause an increase in the lithostatic 
and hydrostatic pressures With depth. Other initial 
conditions cculd be incorrforated into the mathematics if 
required at this stage. Since the compression é =-V 2 E&, 
then at t=0: 
Oe Get neato 
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With the zero initial pressure condition this requires: 
ST Renee 
dz S(t=0) = 0 


The equations simplify to: 
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with p2=k2+S. The solutions must be bounded at infinity so, 


ty the Method of Variations (Butkov 1968, p126) 
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A,B,C can be determined from the boundary conditions. They 


may depend on p,q and s but are independent of z. 


At the free surface, z=0, 6 = @ =0 reduces to: 
XZ yz 
aS Se Ca, -3.15 
The surface oe must equal the weight of the water and is 
represented by the compressive load - T(P,sde0,S)- 
Zl 


aac aae k (u-k)B + Ck (1-a, ) 


The third conditicn involves the volume of flow at the 
bottom of the reservoir or the fluid pressure. Since the 
fluid pressure iS easier to monitor by test wells over a 
large area this was chosen to be the surface boundary 
condition. Water iosses sould not be ignored since leakage 
can be significant. At the Tel-Yeruham reservoir in Israel, 
losses ranged from 10 to 30cm/day over an area of 390,000 
Square metres at the 50m high water mark (Aisenstein et.al. 


(1957)). From (3.10) : 
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From the equation (3.17), (24) 57(3210)) and) (3.14) it is 
possible to compute the stresses, displacements and excess 


pressure at any depth at any time. 


Substituting for E and § into equation (3.10) gives: 
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The terms may be broken into suitable parts by partial 


fractions and the initial values set to zero. 
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The inverse Laplace transforms were all obtained from the 
tables of Roberts and Kaufman (1966). The necessary 
transforms are listed tLelow. Equations (p257, 5) and 
(p224,23) were rederived since they were contradictory in 


the small z limit. 
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evair 2 azt ee hs 
: (Zasctazt lye © mr ect =a) (p257,5) 
(s *+a) eee 5 Pe 
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f (sta) e T Grite (p169',3) 
Using (3.13) it is rcutine to show: 
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From the surface load and pressure at time t, F can be 


calculated by taking the Fourier transform at known times. A 
linear interpolation was used to calculate the wavelength 
contribution at times intermediate between those available. 
Spline interpolaticna was attempted, but proved less 
successful in treating rapid water-level changes. 
Numerically unrealistic results due to overshoots were found 


in the time derivatives in these instances. 


The convolution integration was done using a Gaussian 
Quadrature integration technique (Ralston 1965). The time 
interval was sampied for orders 10,15 and 25 of the Legendre 
Polynomials. Rather than increase the order of integration 
to improve accuracy, it was found that excellent accuracy 
consistent with low cost could be obtained by interval 
integration. Since the derivatives are constant between the 


kncwn depth-time values, each of these intervals was 
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integrated with up to a 10th order polymonial. This has the 
advantage that a small sample interval was used where the 


functions change most rapidly. 


Heaviside loading, or rapid filling is easily 
incorporated into the analysis. If the lake fills at an 


infinite rate to the depth described by I, and is held at 


that pressure for all times greater than t=0 then: 


d 2 — Bs — 

qr r= 6(t) F(t=0) 
t . 

a7 a: 
{& Tr Cte) et (2) dt is PO) EE) 
fe 


In the event of rapid filling the convolution becomes a 
simple multiplication of the [ (just after filling) with 
the function evaluated at time t. 


3.4 Displacement 


The vertical displacement can be found from: 
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This expression may be expanded by partial fractions and the 
inverse Laplace transforms taken from those tabulated 


earlier. 
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3.5 Effective Stress 


The effective stresses are obtained in a similar way . 
The boundary constants A,B and C are substituted into 
equation (3.10), the stress-strain relation for effective 


stress. Expansion by partial fractions and use of the known 


Laplace transforms leads to the required results. It can be 
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The expression for 99 Can be obtained by replacing all 


occurrences cf £ with g in 3.21. 
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aol 1 eee 
Oey [k Erfc(-kt ”) ue ho eee Erfc (At?) ]} -3.23 
zk (a,—1) 
where ag = (2a,n-1) 
pes can be oktained by replacing p with gq in the above 
equation. 


3.6 Vertical Flow Velccity 


The flow velocity used in D'Arcy's Law was defined in 
equation (3.7). It is of interest to know its surface value 
Since it could have been used as one of the surface boundary 
conditions. Substituting into (3.10) the values of A, B, and 


C leads to: 
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where K is the permeability of the rock matrix constituting 


the half-space. 


To illustrate an application of the equations (3.18) to 
(3.24), a two dimensional cross-section of Lake Kariba, 54m 
deep, will be analysed. The two dimensional analysis was 


used only for simplicity in representing the results and is 
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obtained by setting q=0. It should be emphasised that the 
technique is general, and some three dimensional examples 
wiil be given later. To avoid complications arising from a 
time dependent history of filling and emptying Heaviside 
loading wili be used. This may be interpreted as a rapid 
filling to maximum depth. The elastic parameters of the 
drained sample were the same as those used earlier with 
Poisson's Ratio 0.27 and Young's Modulus 0.85 Mbar. The 
other two parameters Q and a were unknown for the rock so 
it was assumed that it behaves as what Biot calls a 
"Saturated Clay*, and which Jaeger and Cook (1976) indicate 
applies to rccks under lakoratory conditions. 
O08 a= 1 

No data was available for the permeability in this 
region so a value of 2 millidarcy with a water viscosity of 
0.01 Poise was chosen and iS consistent with the chosen 
elastic constants. This value corresponds to a sound rock of 
low permeability. Scholz et.al.(1973) has estimated a 
permeability of 2 millidarcy from precursors time scales 
assuming diffusicn through a solid matrix. This permeability 
could change by several crders of magnitude if the material 
is extensively fractured, but at the depths of 1-10 km where 
the seismic energy is thought to be released, these 
fractures are probably not a major factor in bulk fluid 
transfer since they will te closed by the weight of the 
overburden. The permeability of sound rock is most likely to 


describe the diffusion, but it shouid be emphasised that the 
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behaviour of acquifers or major faulting has not been 
included in the model and these may seriously distort the 


diffusion results. 


The pressure at the surface was also unknown. The 
filling of the reservoir changes the water table but at this 
stage the amount of this change over a large area is 
unknown. The literature has not provided a guide to the 
solution of the question and a physical argument must be 
used. If the filling of the reservoir does not alter the 
water table and the surface pressure is zero at all times. 
Examples are given cf the pressure and displacement in 3.1a, 
3.1b. The load in this case is equivalent to a lake being 
filled with the water in the lake separated from the water 
in the rock. This would be the case if the lake bottom was 
sealed or the lake was a biock of ice. This has been 
designated as zero coupling to indicate that the water table 


has not changed. 


Figures 3.1c, 3.1d show the results at various times 
alceng an axis directly below the centre of the iake. Notice 
that the displacement increases at all points with time 
after an initial elastic compression. This consolidation is 
analogous to the behaviour cf a loaded wet sponge. There is 
a certain amount of initial compression and then a_ settling 
as the water flows out of the sponge. The pressure 
variations are particularly interesting. The surface 


pressure has been forced by the boundary conditions to be 
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Figure 3.1: 


For the permeable fluid filled half space 
(a) shows the vertical deflection and _ (b) 
the pressures below a two dimensional load 
applied to the surface 1 year after 
loading. The excess pore pressure at the 
surface is zero at all times. The maxinuno 
deflection and excess pore pressure under 
this lcad are shown as functions of the 
depth below the surface and time in 
figures c, d, e, f£. The half space has a 
Young's Modulus of 0.85 Mbar, Poisson's 
Ratio of 0.27 and a permeability of 2 
mililidarcy. 
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zero, but everywhere else the pressure rises. The loading 
has compressed the fores and the pressure of the water has 
risen to oppose the ccmpression. The smaller pressure 
increases at greater depths is expected since the 
compressive effect of the load decreases with depth. The 
pressure gradients cause flow of the water in a manner that 
will reduce the gradients and this is indicated by the lower 
pressures at later times. The displacements and fressures at 
several depths kelow the surface are plotted as functions of 
time in figures 3.1e, 3.1f. The consolidation of the 


reservoir is clearly seen in figure 3.1le. 


A second possible ccndition that may be imposed on the 
surface pressure is that in which the water pressure rises 
everywhere to equal that at the bottom of the lake. This is 
the expected conditicn for the appication of a water load 
and is designated as ‘full coupling' of the fluids in the 
rock and in the lake. This is indicated by the unit coupling 
on the diagrams. In this case figures 3.2a, 3.2b show there 
is no time variation of either the pressure or settlement. 
The pressure gradients established are such that the flow of 
water in the rock is governed ky the rate of change of 
strain rather than the rate of change of pressure in a 
normal diffusion equation for the material constants chosen. 


This is seen from equation (3.9). 


In this problem the pressure increase was expected to 


be delayed some time tkLehind the pressure increase at the 
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Figure 3.2: 


For the same half space and load as in 
figure 3.1, the excess pore pressure and 
deflection are plotted against depth for 
the surface boundary condition that the 
pressure is the same as the applied 
compressive load. Notice in this case 
there is no time variation. 
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Surface . For this reason it is necessary to examine the 
assumptions made. The consolidation at Kariba was small, and 
as the elastic constants appear to give a reasonable 
agreement between the calculated and cbserved deflection so 
the constants seem reasonable. The three fluid parameters 


chosen for the calculaticn are the permeability, oF , and Q. 


The permeability only alters the time scale, so to 
expand the time effects a low permeability has been chosen. 
Increasing the permeability would only force the fluid 
equilibriun to be established earlier. It should be 
appreciated that if the reck is fractured extensively, the 


fluid may move through the fissures. 


The seccnd parameter o% relates to the rock saturation. 
It is reasonable to assume that rocks are relatively 
saturated so that os lies between 0.9 and 1. A larger value 
of om decreases the time required fcr equilibrium, however, 
if ™ is too smail the pressure may increase in the manner 
shown by figure 3.1. This is callled the Mandel-Cryer effect 
and the ancmalous pressures may be larger than thethe 
applied pressure ( see Schiffman et. al. (1969)). For large 
Q, noticeable time lags can be seen for o& of 0.9. The 
consolidation in this case is very small and is iliustrated 


in figure 3.3. 


The other parameter that may be varied is Q. This nas 


been specified in terms of the initial (undrained) and final 
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Figure 3.3: 


The displacement, excess pressure, and 
stresses directly below the two 
dimensional lake section at several times. 
The elastic constants were the same, but Q 
and alpha have been reduced to illustrate 
that time delays by diffusion are 
possible. 
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(drained) ccmpressibilities. In the laboratory results 
quoted Ey Jaeger and Cook, a rock with high compressive 
strength will consolidate by the movement of water. The 
final compressive state is almost entirely a result of the 
water volume change. For this case Q is infinite with the 
final equilibrium being established instantaneously. Figure 
3.4 illustrates that lowering the Q value leads to diffusion 
and a time lag. In the geometizxy of figure 3.2 detectable 
time lags at 0.1 years requires that Q be about 10 times 
larger than Young's Modulus. This implies that the initial 
to final compressibility be in the ratio of 1/10. The 
surface consolidation would still be smaller than the 


accuracy of relevelling data. 


The literature on Q and om for _ rocks is limited. 
Hendron (1969) presents the results of compression tests on 
drained and undrained black shales. From the Young's Modulus 
found, it is possible to ccnclude that Q is of the same 
order of magnitude as_ the drained Young's MOdulus. 
Computation with a Q five times larger than Young's Modulus 


leads to the results shown in figure 3.4. 


With the variation of the elastic and fluid parameters 
of the rock a large range of results can be obtained. A lot 
of work should be done to determine these parameters before 
the numeric results are emphasised. The caiculations made 
here will concentrate on the behaviour, rather than the 


magnitude, of the rock stresses. 
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Figure 3.43 


The displacement, excess pressure, and 
stresses belcw the two dimensional lake in 
a fully saturated rock. The elastic 
conditions have remained unchanged, but Q 
has been reduced to illustrate that 
pressure increases with time are possible 
with a different choice of fluid 
parameters. 
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Another approach to the problem is to assume that the 
surface boundary conditions are poorly chosen. If the 
surface is heavily fracture, inhomogeneous, and geologically 
variable, then the protlem should be initialised at some 
depth just below the surface. At this point the vertical 
load is still similar to the applied ioad but the water 
pressure is assumed smaller than the full head. The boundary 
pressure would probably extend beyond the edges of the lake, 
but for simplicity is assumed to be of the same shape but 
smaller by a factor of €.25 than the applied load. This is 
assume to be 0.25 coupling, to imply that the load lies 


somewhere between the blicck of ice and the fully commected 


water solution. In the absence of odther vaiues Q and were 
chosen to be , and 1 respectively. The theoretical [f now 
becomes Ee x Z (2a,n-1) es 
T = = — a 
(a,T + aP) ( (-a,) +9 0. 250) 0 


where a =0 for Q= © a =1 
Figure 36:5 summarizes the pressure and displacement 
variation with 0.25 coupling. The settling is smaller than 


with zero coupling. 


The equations (3.18) to (3.24) were analysed using time 
convolution. A typical loading history was modelled and is 
shown in figure 3.6. The water depth relates to the lake 


cross-section shown in the two dimensional sections of 
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Figure 3.5: 


The excess pore pressure and displacement 
are plotted for a permeable half space 
with Young's Modulus 0.85 Mbar and 
Poisson's Ratio 0.27 and permeability of 2 
willidarcy for the situation where the 
water table has increased by 0.25 of the 
depth of the lake. Figures c, d, e, f£ show 
the variations of the maximum deflection 
and excess pressure as functions of depth 
and time. 
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Figure 3.6: 


Plot of a hypothetical depth filled as a 
percentage of the maximun depth against 
time for the two dimensional lake section 
shown earlier. The downdraw starts at 4 
years at 90% full, drops to 50% at 4.5 
years, and from 5.0 years rises to be 100% 
full at 5.5 years. This filling history 
will be used with other lake sections 
later. 
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PVGUneSeaw let. 2. This Same filling histcry , although to 
different depths, will be used in the other examples at 
Kariba, Nurek, and Oroville. This variation of the water 
depth was meant to illustrate a typical downdraw situation 
where the lake is drained for a short period for dam 
inspection and then rapidly refilled. There was a similar 
downdraw just before the seismic activity began at Lake 


Oroville ( Marlette 1975). 


To examine the effective stress just before, during and 
after the downdraw, plots showing variations with time and 


depth were made for the loading history in figure 3.6. 


The figure 3.7 shows the cross sections of pressure and 
displacement at 3.75, 4.75 and 6.5 years at Kariba. The 
pressure at 3.75 years is less than that determined in the 
earlier results shown in figure 3.5. This is a result of the 
Slow loading of the lake which allows for the fluid pressure 
to dissipate as filling progresses. Instability is difficult 
to define, and relative Mchr circle positions were used to 
decide when an area is ieast stable. If at any time a Mohr 
circle lies to the left of its position at another time, or 
is larger and hence closer to a failure envelope, 
instability is suggested. The stresses produced by the lake 
loading are relatively small, only a few bars, and must of 
course be added to pre-existing tectonic stresses of up to 
several hundred bars before the Mohr circle is meaningful. 
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Figure 3.7: 


Cross-sections below the reservoir at 
3.75, 4.75 and 6.5 years showing the 
excess pressure and consolidation. fhe 
parameters used are shown on each diagran 
and they were all computed for 0.25 
coupling. 
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at Kariba where dip-slip faulting occurred a will equal 
the weight of the rock and ae must be less than ee A 
review of the state of the stress in the earth was published 
by Haimson (1975). It appears that the horizontal stresses 
may exceed the vertical at shallow depths (1 km) but are 
jess than lithcstatic below 1.5 km. The ratio of horizontal 
to vertical stress being eventually about 0.6. These 
conclusions were broadly applicable to southern Africa and 
the United States. It of course does not apply to thrust 
environments where the largest stress must be horizontal. 
There is little reason to expect these cbservations, which 
are averaged over large areas, to apply to a particular 
place. At each dam site under investigation an attempt to 


determine the local tectcnic stresses should be made. 


Figure 3.8 shows the variations of stress, pressure, 
and displacement as functions of time at: a point centred 
under the load a depths of 3.,6., and 12. km. The sum of the 
stresses indicates how the Mohr circle moves to the left and 
right and the difference is a guide to its radius change. 
The Mohr circle position and radius change with time and 
will be examined in a later section. Figure 3.9 shows the 
changes of variables with depth at times just before (3.75 
years), during (4.75 years), and after the lake downdraw 
(6.5 years). The rebound is more than expected but reloading 
is sufficient to consolidate the material more than that 


observed before the downdraw. A different loading history is 
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Figure 3.8: 


The variaticn of the vertical 
displacement, pressure, and stress in the 
Z and xX direction, and their sum and 
difference are plotted at 3,6, and 12 km. 
The displacement (a) illustrates how the 


unloading relieves a lot of the 
consolidation and this is reflected in the 
pressure low in (b). The pressure 


increases initially but as the rate of 
loading decreases the pressure falis as 
the fluid diffuses away. The changes in 
the effective stresses (negative in 
compression) were staller than 
anticipated. Notice the change in the sun 
of the stresses occurs later than the 
Major change in the difference. 
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Figure 3.9: 


Changes in displacement, excess pressure 
and effective stresses are plotted for 
different times as functions of depth. The 
values chosen for the constants agree with 
those used in figure 3.7. The loading 
history was shown in Figure 3.6. 
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used for the Figure 3.10. The change of rate of filling has 
little effect on the centre of the Mohr circle, but the 


radius follows the level changes. 


Three dimensional sections can be analysed and an 
example is shown in Figure 3.11. This shows the pressure, 
displacement, and effective stress at 5.5 years after 
filling using the filling history of figure 3.4. It is felt 
that the apparently noisy results are real since similar 
results were cbserved at higher orders of integration in two 
dimensions. The stresses are difficult to interpret in three 
dimensions, so an analysis of the physical behaviour will be 
done for two dimensional cross-sections only. The figure 
3.11 indicates that the same results will generally be 
applicable. A similar analysis was repeated for lakes Nurek 
and Oroville. The results are shown in figures 3.12,3.13 . 
In each of these cases the filling history was the same, but 
the depth of fill varied in each case. The bathymetries are 
given with the lake depths being found by the technique of 


triangular interplation cutlined in section 5.1.4. 


3.9 Layered Porous Flastic Half-Space 


To allow for the variation of parameters with depth two 
approaches are fossible. A functional dependence, perhaps 
exponential, may be introduced and the derivations found 
analytically. Substitutions must be nade into the 


expressions already derived and similar processes followed. 
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Figure 3.10: 


Displacement, excess pressure, and 
effective stress for the loading rate 
shown (a). The changes in filling rate are 
Only slightly seen in the movement of the 
centre of the Mohr circle. 
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Figure 3.11: Three dimensional analysis below Lake 
Kariba at 5.5 years. 
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Figure 3.12: Three dimensional analysis below Lake 
Nurek at 5.5 years. 
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Figure 3.13: Three dimensional analysis below Lake 
Oroville at 5.5 years. The excess. pore 
pressure, vertical deflection, and 
effective stress are shown. The loading 
history of figure 3.6 was used for the 
calculation. 
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The algebra will become even more confused than it was in 
the work shown above. Ancther approach models the halfspace 
as a series cf layers, each with different elastic and fluid 
properties. This is often a gocd approximation to 
sedimentary basins, and if there are a sufficient number of 


thin layers reasonable parameter variations may be modelled. 


The treatment is analogous to that presented earlier 
for the elastic case with each interface having six boundary 
conditions. The continuity of displacements u, v, w reduce 
te two conditions which requires the interfaces to stay 
together. Continuity of the three stresses a again reduces 
to two conditions which prevents acceleration of the 
interface. Two cther constraints are required. The pressure 
of the water must be continuous across an interface and the 
amount of water flowing out the bottom of a layer, ae must 
equal the amount flowing in the top of the layer below it. 
The time scaling introduced earlier depended on permeability 
and cannot be used here with each layer having a different 


permeability. 


The equations (3.11) and (3.12) apply at any point. 
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equation to: 


2 d a 
Car oe <5) (-k* + Seti =- 2a,s <s 
(oA az ? 

ae. 
(eee ee 
az 
where = K2 a s 
Cc 


Solving the differential equations leads to the expanded 


form: 
a teks ames += = = 
12, Se Nae =- JN ES . + Be na + B ee + a,c ze kz = a,c eS: 
“n~ ate = = 
Sy ay a ara 


: + - + ~ + - 
The six unknown constants A, A,B,B,C andcC are 


determined from the continuity conditions at each interface. 


Substituting into equation (3.10) leads to the matrix 


equation: 
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and 0 is the 3x3 zero matrix. 


The stresses are determined from the matrix equation: 
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Continuity at the i interface at depth z below the surface, 
separating the i and i-1 layers reguires: 
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JE 1 


or 


[al,_, = Petz iJastmlsc,, (CMa; boustal ay a beizy)ilals 


—3).27 
The notation cf figure 2.9 has been used here. Below the 


lowest interface (N) lies a halfspace and the boundedness of 
values requires a B and eS to be zero. The top surface 


has specified values of excess pressure, 5 , the normal load 


co and oe =0. The three remaining unknown constants can be 


evaluated from: 
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fhe equations (3.27), (3.28), (3.29) are sufficient to 


calculate the A,B, ey in all the layers (i=1,2,...N) 


+ 
a: =a 
and (3.25) and (3.26) suffice to compute the stresses and 


displacements everywhere. 


The inverse Laplace transforms are difficult to do in 
this layered model. The inverse matrices puqeas can be found 
é 


analytically in terms of u and when post multiplied by 


Cay, 2 in the next layer leads to functions for which the 
’ 

inverse Laplace transforms are known. After several 
attempts, effcrts to find Simplified inverses were 


discontinued. There are algekra routines (REDUCE2, ALTRAN, 
FO RMAC) available which can do this type of tedious 
Manipulation but they were not used. It is sufficient to 
conclude that [A] in an upper layer may be connected to {Aj 
in other layers by a series cf convolution integrals. The 
convolutions become derivatives when a bottom iayer is 
connected with a top layer. This series of time convolutions 
and derivatives indicate the results at any time and place 


depend on the behaviour at ali other layers and times. 


The numeric instability inherent in the required 
computation indicated the results would always be suspect. 
The analysis is adaptatle, but other methods are probably 
more useful for this type of problem. The finite element 


technique will be described in the next chapter. 
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Chapter 


i= 


Extension and Improvements of the 


Work 


The theoretical descripticn given in Chapter 3 is only 
a first attemgt to describe the behaviour of water filled 
material. Various extensions can be made to the theory to 
allow for changes in the various parameters which depend on 
stress. Assumptions have also been made which, though 


unnecessary, simplify the model representation. 


The stresses which are computed are only the anomalous 
stresses applied to the ‘neutral* background with the 
initial stresses, hyroicgic conditions and the gravitational 
forces being ignored. These are easily substituted into the 
equations but were omitted to avoid confusing the output 
results. In any applicaticn requiring the absolute stress 


levels these cf course shculd be considered. 


The simple, or conventional effective stress law was 


defined in section 1.3 with a=1 being used in the equation 


° = Oo ; “ym eae 
OC affective applied 


where P is the fluid pressure. For large changes of stress 


a is fcund to be pressure dependent. Terzaghi argued that 
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(= ne where 1 is the porosity, however, his experimental 
results Showed that 4% was close to one. Handin et.al. 
(1963) extended these experiments to different pressures and 
rock types and showed © is usually one except for very low 
porosities. Brace and Martin (1968) try to explain this to 
be a result of the applied strain rate. Brace and Martin 
also refort a pressure drop (dilation) just before failure. 
Hubbert and Rubey (1959) attempted tc preve that o=1, but 
Nur and Byeriee (1571) show this argument is not valid for 
low porosities. Skempton (1961) and Geertsma (1971) proposed 
a model which was verified by Nur and Eyerlee (1971) in 
which 

ao = 1- B/B, ae 
where 86 is the bulk modulus of the porous rock and Bs is 
the bulk mcdulus of the rock grains. Usually B/B | is small, 
but 8 is strongly dependent on pressure. Empirically they 
suggest a linear dependence: 

bi =1 Bid Pee gee) -4.2 

where Pe is the water pressure and P the confining 
pressure. Knutson and Eohor (1963) show the Variation of 
bulk modulus for several rock types over a wide pressure 


range. Changes of more than 100% are observed for Weber 


Sandstone with a pressure change of 1000 bars. 


Moreland (1972) apprcached the problem of defining the 
pressure and stress as the average Over several grains. This 


work on the ‘Theory of Interacting Continua (TINC) derived 
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the form of (4.1) by a first order simplification. Garg 
and Nur (1973) compared each of the models and concluded 
that while the conventional stress law (a=1) overcorrects 
for the pore pressure effect, the TINC and (4.1) 
underestimates the effect. They suggest the underestimation 
may be a result of strain rate used. They also show that the 
conventional stress law best describes the material strength 


but is the poorest in describing deformation. 


This study is mainly ccncerned with fracturing and the 
conventional stress law (a=1) was used. If the application 
is directed more towards deformation the equations 4.1 and 


4.2 should be used. 


The existence of fractures reflects the water 
transferring capability of rocks. If an earthquake occurs in 
a certain region it will introduce new fractures and may 
change the drainage pattern of the area. Water may also be 
transmitted by diffusion through the rock, so _ the 


permeability and fracture frequency must be considered. 


Little information is available about effective 
permeabilities on the large scales of 10's of kilometers, 
however both aperture aspect and permeabilities are stress 
dependent. Snow (1970) describes fracture frequency and 


aperture width at depths of about 25-100m. In most cases, 
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both the frequency and permeability (which depends on the 
aperture) decrease with depth. Fatt (1953) has shown that 
the permeability decreases with increasing pressure in 
hydraulically tested samples. For rocks taken from oil-gas 
reservoirs, the permeability was observed to decrease by as 
much as 30% upon application of 300 bar confining pressures. 
The rocks he tested had high permeabilities initially: 
between 100 and 250 millidarcy. Morgenstern and Guther 
(1974) ainterpret experimental results on dry samples 
Slightly differently. The permeability was certainly found 
to decrease with fressure, it decreased rapidly initially 
and less rapidly at higher stresses. In most cases the 
permeability of sandstone was reduced by less than 50 per 
cent over substantial pressure ranges. Coais, granites and 
gneisses had larger reductions attributable to changes in 
the secondary structure of the materiais, Guther (1972) 


found that the function 


embraces much of the experimental data. He defines the 
average principal stress 
ik 
_i 2) (1a) (to, hoe). | 
J = 5lo,(1+2R,) (1+v) (Ao, +40, 
Where. 0.) co. , Lo -, Ac are the horizontal and vertical 
Vv h Vv h 


effective stresses and the changes in the stresses 


respectively. R is the ratio of horizontal to vertical 
fe) 


effective stress. ) is Poisson's Ratio. "A* is a constant 


adjusted to match the permeability at zero effective stress. 
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The values of T and N are found from experiments. For rocks 
Showing a strong dependence of permeability with stress, the 
best approxiwaticn for the samples tested was N=3.0 T=20 psi 
(2.8 bars). Weak dependence gave N=2.0, T=40psi. These 
values will depend on the rocks used but indicate the 
typical magnitudes. It shculd be mentioned the fit of this 


function was shewn only in the 1 to 25 bar range. 


Gangi (1976) has derived a permeability-pressure law 
using arguments of Hertz compression of spherical granuiar 
material. The permeability variation with pressure is given 
b 

Y (P+P. ) ees 
=| ) 
I 


Oo 


Coe = Lyssa - om 


where Ko is the initial permeability of the loose grain 
packing, Co is a packing constant (of the order of 4), a is 
the “eguivalent pressure due to the cementation and 
deformation of the grains since there is some deformation 


even at zero pressure. as is the effective elastic modulus 


of the grains defined by 


ze 2 
et lea) 


where 8 is the grain*s bulk mcdulus and vy is Poisson's 
Ss 


Ratio. 

Gangi (1976) also derived the permeability-pressure 
relation for a "bed of nails" model. In this case two flat 
layers are separated by small, column like, asperities. The 


functional dependence iS K(p) = K (1 - (P/P,)™)° 
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where KS is the zero pressure permeability, ery Be 1 ey 
effective modulus of the asperities defined by 

Py = 5 AL/A ie ts 
and eee is the fraction of the fracture face covered with 
the aperities (generally this is small). ‘m* is a constant 
(O<m<1) which characterises the distribution function of the 
asperity lengths. Gangi was able to obtain a good fit to the 


experimentally determined fermeability of a 200 millidarcy 


sample between 0 and 700 kars. 


In geologically realistic probiems the earth is rarely 
homogeneous Or isotropic. There are always rock type 
changes, fracturing and large variation of parameters. fhe 
largest single disadvantage to the thecry outlined is its 
inability to handle lateral variations of parameters - and 
this includes inclined faulting. One way to solve these 
problems uses the technique known as ‘Finite Element 
Analysis'. In this, the space is divided into elements and 
the behaviour within the element is found by a variational 
analysis with the value at the nodes either being specified 
or computed. The determination of unknown nodal values is 
aoe by solving a matrix equation. Within each element 
interpolation is done by folynomial expansion, usually of 


orders one to four. 
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Finite element analysis is not often used simply 
because a three dimensional froblem requires a prohibitively 
large matrix inversion. Problems with axial symmetry or with 
two dimensions are quite solveable using this method. The 
time dependence is another difficulty since linear time 


stepping introduces numeric noise. 


For completeness a description will be given of the 
development necessary for finite element analysis of Biot's 
Consolidation Equations. The coupling of transform 
techniques with F.E.M. for solutions of the equations near 
the surface shows great promise. Near surface fracturing, 
rock changes and topography couid be soived by the F.E.M. 
and below this the analysis could be done by transforms 


since the prceperty changes are more gradual. 


The literature on the F.E.M. technique being used in 
consolidation problems is indeed extensive. A selection of 
the literature is given in the Bibliography with the works 
of Hwang, irobe, Valliappan and Zienkiewicz being 
noteworthy. The notation and description by Irobe et.al. 
(1974) will be followed since it describes unsaturated 


soils. 


4.3.1 Variational Equation 
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The basic equations are: 


(a) the equilibriun equation 
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where © is the stress and P, the density of the saturated 
1) 
sors F. are body fcrce ccmpcnents with gravity being the 


Significant force. 


(b) the constitutive eguation 


AR eet ae 
wees Ci 5ke “ke 


where the Ba are displacement components of the solid. 


Ghaboussi and Wilson (1973) show that the introduction of 
two elastic matrices will describe a compressible fluid 
model. The 2g is a second order symmetric tensor and 

are the 21 elastic moduli of the solid. 'P* is the excess 
pore pressure. 


(c) Strain displacement equation 


iu 
= ssl. 4 ap Wha x 
5) 64s air 
(d) the generalised D*Arcy's equation 


3 = a = = KK. « eet: F -4.3 - 
MS OY ler ee Pobe! 4.3 


where Ve is the component of fluid flcw volume, ae is the 


permeability tensor, p, the mass density of the water. U 


is the compenent of fluid displacement and u, the 
displacement component of the solid. The inflow volume ile 
rel 


is composed of air void, volume change and solid 
deformation. If the Bict notation is again used, 1/Q0 is a 
measure of the amount of water that can ke pressed into the 
solid while the volume is kept constant. Then: 

ce Vig & eee a Rs array 
Two invariants must be used in a variational solution to 


this problem. The first invariant is defined: 
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and it may be referred to as the elastic consolidation 


potential. If the initial conditions are all zero: 


= = Iu). ee et: = 
u; 0 U; ete Js 0 


then (4.4) may be rearranged: 
1 -4.5 
== : pn 3t Uni : 

a AE WA al Naa cl a ch 
Le Q plays the role of the dissipation function, then its 
variation is given by 

= ay Vv. 6V.dv 

62 Hie ee i 

where Noe is the flow resistivity tensor equal to the 
inverse of the fermeakility tensor. The variation equation 


must also satisfy the boundary condition and the total 


variation for the problem becomes 


6¥ + 6Q = fea Fe Su, dv - i ie OV. dv 


~ 6V n.)ds =H 0 
+ fey Sts SES n;) 
where 7T. and P_ are the surface tractions and pressure 
a Oo 


acting on the surface with outward normal ° 


If the space below the load is divided into elements 
and the value at any point is described in terms of a 
polynomial expansion cf the nodal values, then the 
displacements and inflow volumes may be expressed in the 


Matrix form. 
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Capes (een 2185 


{v} = [nN] {v}> -4.8 
where the notation { } is used for column matrices, < > for 
row matrices and [ ] for rectangular matrices. The [NJ] are 
the shape functions for the elements chosen and { }* implies 


that the array contains values of nodal foints. 


The strain matrix can be obtained from the derivative 
of (4.7) 
{e} = (B] {0} * -4.9 

where ([B] is obtained from {N] and the equation (4.3). The 

. pressure is determined from (4.5), (4.7) and (4.9). 

Pp = -9(<1>[B,]{v}" + <A>[B]{u} ) 
where <1> = <1,1..,1> and < > = <a,1,-..%3;>. The [B, } 
matrix is required for contraction of the indices in the 
{ }*. The six components of equation (4.6) may be written 
in matrix form: 


SY fsu}T[B]2[p] {ul av 


* 
+0 f cor*™ (te? <a> De) a)” + [B]7{A}<1>[B,]{v} ) av 


* iy * 
£0 few"? (a, a3<0 [elo ol beall Clie le (iv) dv 
-4.10 
sa = [ tév} “[N] "(pI IN] {viv San 


where [D] and [D ] are the elastic and flow resistivity 
V 


matrices for the elements, and 
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Tig}* 


fe 6u. ds = fu} {F} =A, 14 
ak Alb 
fe Su ae = tee) Pai econ 
ES al O 
-4.15 


where (PJ is the column vector of the body force ae, ’ 


{F}* is the cclumn vector of the equivalent nodal point 
loads, [({H]) is the matrix derived by the integral of the 
product of the interpolation functions of inflow volume and 


pore pressure along the boundary, and a ie are pore 


pressures at the ncdal points. 


Since the variations of 4 and V are arbitrary then 
substitution of (4.10) to (4.15) into (4.6) leads to the two 


equations describing the equilibrium equation and the flow 


equation a * * * 
Phe Ge ([K]+[L]) fu} +[L, ]{vi - 1Fj} -4.16 


* * * e * 
Se eS [L,]ful = [L,]{v} + [M]{v} - 1F 5} -4.17 
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It is assumed by Irobe et.al. that the v changes linearly 
over small time intervals A so that { Vv 3}* evaluated at 
time t may be expressed as: 

* 2 es 2 
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2 x } 
= mn wire om CG t=-/A 


Equations (2.16), (2.17) may be combined in the matrix form 


r]> be i ig bs ; a igs ae U 
=c : F, . Gal ee Kia Koo V : 
where [k,,] Set speed 
[kK,5] = [4] 
[k,,] = [t,] 
ifr genie cei ntl 


This equation can be soived by time stepping as it has a 
recurring form. This time stepping is numerically poor and 
may diverge. Other technigues using time convolutions (see 
Hwang et.al. 1972) avoid some of these difficulties. The 
earlier analysis showed that the time dependence usuaily 
involves error functions, however published results using 
time stepping agree well with analytic solutions in two 


dimensions. 
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In three dimensional problems several element shapes 
are possible. Zienkiewicz (1971) describes these elements 
and indicates that the choice will depend eventually on the 
problem. Since there is nc symmetry in the lake problems the 
choice would be limited to the tetrahedron and the brick 
element. The number cf ncdes chosen will also depend on the 
order of polynomial the interpolation will use. In general, 
polyncmial order wili ke one less than the number of nodes 
along the element edge. Examples of elements are given in 


figures 4.1 and 4.2. 


wothe tetrahedral element shape functions are best 
described in terms cf volume coordinates. These coordinates 
have values between 0 and 1 and are equivalent to the 
normalised perpendicular height above a particular side. 
There are of course four coordinates for the four corners 


and heights. In figure 4.1d the coordinate L; is defined by: 


_ Volume (P234) 


GS volume (las ap 


For the cubic tetrahedron (4.1c) the shape function for each 
of the node types is: 
Corner Nodes 


Midside Nodes 


Mid-face Nodes 
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figure 4.1: 


The tetrahedral finite element family with 
the required nedes for (a) linear, (b) 
guadratic, and (c) cubic interpolation. 
The volume coordinate, with values between 
0 and 141s defined as the ratio of the 
volume (P234) to volume (1234). 
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figure 4.2: 


Nodal points for the ‘brick elements. The 
serendipity elements (a), and the 
Lagrangian elements (b) are shown for 
linear,guadratic, and cubic interpolation. 
Notice the large number of internal nodes 
required for the Lagrangian elements. 
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The integraticn fcrmula required in the integrations of 
(4.18) is 


aye Dc dd a!lbici@d! 
L [ir Jie = ee ee eee ee 
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element 


The matrix fcrm of equation (4.7) can be written immediately 
with this knowledge of the shape functions N . Since all the 


nodes there are three displacement components specified: 


d 36 
u = [IN,, IN,, «-- rIN5] ay = [N]{u} 
V 
\ rey 
WwW 
eh 
lee 
5 
ea 
where I = il 0 0 WwW. 
0 1 d 
0 0 


In the case of rectangular prisms two shape functions 
are usually chosen. The ccordinate axis, shown in figure 
4.2a, is centred in the middle of the brick with the bricks 
edges having coordinate values of +1. The cubic 
"serendipity" elements require nodes only along the edges 
and the polynomial expansion is incomplete. This is still 
satisfactory since the F.E.M. usually only requires 


continuity of the boundary values and the first derivatives. 


The shape functions are derived in terms of the 


coordinate axis given in figure 4.2a. As an example a cubic 
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element has 32 nodes and the shape functions are given by: 


Corner Nodes 


ii pee ea 
Ne Sean Clipe) Chama il he +n +6) ] 


Mid-side Nodes 
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These forms were introduced to allow for + signs in the 
general formula. The other shape functions useful for brick 
elements involve products of Lagrangian polynomials which 


can be written 


This is clearly equal to one at the point x . and is zero at 
all other iy The shape functions in three dimensions are 
products of those in one dimension: 

Seen = L. (é) Bh Ly (v) 
The Lagrangian family is unlimited and is easy to generate, 
however it is limited by the large number of internal nodes. 


The higher order polynomials also exhibit poor curve fitting 


properties. 


The Lagrangian and Serendipity elements in two and 
three dimensions are shown in figure 4.2. The quadratic 


Shape functions in two dimensions are shown in figure 4.3. 


ns fi ( : 
tema” ye et pein, 


i , ©¢ ‘ F + i? 
pare? ha +t? te oiy aa a 
er Sees | ~ a 
= * oY * ; ; = 3 
a" Pat i = =) a sila g 


ai? AL sopke £ 20% welin oF banuboxsat 91s" nee 
jnteg sot fvJony apoi7oast saad zedto dT Sivaso? “iwseae 
jckiv @£eie0uyhog osipasaped Te etouborg evlovak s2ae0io 


a. 


asttiiv od aBD 
np A 
is 


. eer “a yeah i ct-x) ar i 
ra duie . i ie St 
ey oe a ROA c *y i $ r 
i< 
fs onss ab bas, * aPt0e “ ts spo ua Larps gixeelo a mat 


sie <bdiddaaks. sath? oe saoitonud aqade odT 2 sodso Ltn 


eg kenpas h ono td anaes to eeovboxq . 


oe 0 
a Ms eg | 
ge iz ato ale pes 


a alah ods dees 2i Sue hepeakany ea ean coigeowsl » 


picendiiiaievaie shed-onmpen, ah 


. Setieieay 6 
ir ee a 


= 


= ahesiee | Lenoiapowtd ows a: ‘su eiupl2 
Y ote), YJigibnuatee sit x01 ecottonnt sha led 
eigowsls {4} asipasgp 


2 


figure 4.3: fhe two dimensional quartic shape 
functions for the serendipity (a), and 
Lagrangian (k) elements. 
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4:4 An Introduction to Statistical Prediction 


A large amount of statistical analysis is done in an 
attempt to predict earthquakes and determine seismic risk by 
examining the seismicity of a region. This approach is of 
particular interest in areas such as Kremasta (see figure 
1.3) where it is difficult to determine whether a certain 
earthquake is related to the filling of a reservoir, or is 
part of the natural sequence. Oroville (see figure 1.7) is a 
case where there may have been a change in the pattern of 
earthquakes. The event in 1975 occurred in an apparently 
quiet area, and the decisicn should be made whether this 
fell within the scatter cf the normal tectonic events or 
whether it reflects a change in the seismicity . The answer 
to this question is relevant to other dams in the area. If 
there had been an induced change in the seismicity with a 
delay after filling, then it may happen again. If the 
Oroville earthquake was of natural origin, then other dams 
will have different behaviour depending on their proximity 


to tectonic lineaments. 


In an attempt to answer questions of this type a review 
of the literature was made. At this stage no attempt has 
been made to apply the techniques described below to real 
lake situations and their inclusion is fcr completeness and 


as a guide to future work. 


The usual work on seismicity is described in a paper by 
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McGuire (1976) and involves the determination of seismic 
risk and maximum intensity for state building codes. The 
area under study is divided into regions each with similar 
seismic and geolcgical features. It is necessary to 
determine the magnitude of an earthquake which may occur 
with ¢ specified probability. The events are assumed to be 
Poisson distributed when aftershocks are removed from the 
sequence. The size of successive events are assumed 
independent and the number of events are exponentially 
related to the intensity. McGuire (1976) gives a description 
of the aigebra required, and shows how the uncertainty of 


the parameters affects the statistics. 


Possibly a more suitable statistical approach is that 
proposed by Knopoff, Vere-Jones , Keilis-Borok and others 
using stochastic processes which are reviewed by Knopoff 
(1971). In this 1971 review of the Markov process, a model 
is suggested in which the energy of a certain area is 
increased at a constant rate and released by sudden drops so 
that an energy kalance is achieved. A typical energy buildup 
and release resembles an irregular sawtocth pattern. The 
probabilities can be specified so that the larger the energy 
at a certain time, the more likely there is of a large 
earthquake cccurring. Sfatial as well as temporal changes 
should be included in this model for it to be applicable in 


a particular region. 


A third methcd of predicting earthquakes has recently 
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been suggested by Gasperini et.al. (1976). From the Italian 
data they used they were able to predict that a ‘large? 
earthquake would occur within six years of a swarm in the 
same region. The definiticn of a swarm is based on the 
number of events in one year and a certain area, exceeding 
the average of the data available up to that time. They 
predicted 14 cf 20 earthquakes using this method. The events 
they missed were the smallest cf their defined mainshocks, 
and mostly cccurred befcre 1920. Two events which followed 
On successive years after other mainshocks were also missed. 
They do, however, ccmpare this technique with the model in 
which the swarm and after-shcck occur independently. The six 
year period would vary between regions but is too long to be 


of use aS a predictive tocl. 


Each of these three technigues suffers from similar 
problems which limit their use in the prediction of 
earthquakes at dams. In each, the prediction is based on 
past experience - they are ne processes. The loading 
of the reservoir changes the stress pattern, albeit by a 
small amount, but enough to question the applicability of 
the statistics after filling. The second difficulty is in 
the determining of two dimensional probability functions. 
The probability function must take account of the fact that 
most earthquakes occur near pre-existing faults, so seismic 
lineaments must be given a high weighting. Allowance must 


also be made for a normal spread in the efpicentres, and the 
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uncertainty in determining the hypocentre. However the major 
limitation of these statistics is their lack of resolution. 
It is felt that the use of statistics will never be 
sufficient to decide whether an earthquake will occur with 
the resoluticn cf 10"s of kilometers , and less than a year 
necessary for design considerations or cause-effect 


considerations. 


Statistical research is an excellent field of endeavour 
Since the present techniques seem unable to solve the 
statistical questions raised at dam sites. This work would 
also te relevant to nuclear-site studies where, thus far, 
most of the work has concentrated cn the detection of 


microevents and seismic lineaments. 
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2:1 Conversion of Bathymetry to Regular Grid 


A considerable amount of effort was required to reduce 
a contour map of a lake bottom into a form useable by the 
transform prcgrammes. Three techniques can be used to 


convert randomly distributed data to a regular grid. 


Dr. and Mrs. Gough reduced the contour map of Kariba by 
dividing the map into over 1000 squares. The average 
pressure within the square was found by estimating what 
fraction lay ketween the existing contour levels. fhe 
benefit of this procedure is that good control may be had at 
ridges and saddie points and areas usually subject to 
different interpretations. The process is obviously very 


time consuming and a alternate method was sought. 


5.1.2 Digitising 


In an attempt to increase the speed of the data 
reduction a new approach was used. The data was digitised at 


available contour levels and interpolation programmes used 


han * . : eS 
~~ <a * 


Ale 


be ania 


: ig ups sae 
piistebizacs id 


edt qd eldseal sat 5 oftt 03908) vfs. 0 qe aw0s009 s 
ot beep o¢ ASS eonptaape?> own +2ommerpo3g axoteasad 


enabsy of Beaziupes ase s10%te te | 


2 


sales = : 


Jbi36 ablupss # of stab I - 


ae oe tale 
eet) : 


yd sdigsa Yo gee syosue> Sa? oeoubad sevod -238 boa .20 
epszevs oft .egisdpa. 00h Zewp eras ga ent A ss eign 
dsdy poittentgao ya Bbowod 2st ube edt wiagiw oweeesq F 
eft .elevol wtiso> yaldetep of% aepeged yet aoisosat 
ts Ssd o¢ yew Lozina> boon Sedd ef o1be9039 @idt to sk toad 
ot #oerdee yllewan ssems bas aratoq eibbes bas 2enbi3 
yuev yfevoivde =i a2agoorg sd -2qviseteaysetat saeze22ib | ‘ 
_?epeGe eae bodtew sieuzesis ao base onisveaoo outt | 


| ; = 
stsb sdt 20 baeqe eds senesoul Of #qeerte a8 OF as 


Vie 


46 be@kthokd eau €206 907 -i9e0 exe dpsvdggs wea 6 sotsoubex 
basu cedesryom avisoloqiesar be ajevak =e of pLiw 


199 


to find the value at a given point. at Oroville, seven 
levels each 100 feet apart were digitised while at Nurek, 6 
contours were used with 50 o ena wren Since that contour 
map was less detailed. fhe Oroville digitising rate was 0.1 
inches, however, for reasons outlined later, some sections 
were redone at 0.05 inches when greater detail was desired. 
The number of points (x,y,z) produced by the digitising was 
very large - over 16,000 at Oroville. The large scale Nurek 
map was digitised with a separation of 0.05 inches 
everywhere, but as it was a smaller map, the data points 


numbered only about 8000. 


Several difficulties were experienced in the digitising 
which, although easy to correct, are inconvenient. The 
digitising is usuaily done in several sessions and great 
care should be taken in resetting the origin. The error was 
usually one digitising unit, (0.05 or 0.1 inches) since the 
nearest unit is recorded by the machine. Plotting the 
results as an overlay usuaily makes the error obvious and 
corrections easy. The digitised data is usually stored in 
files, 96 for Oroville, and difficulties arose in 
identifying which file belonged to which contour and where 
on the map it belonged. Eventuaily each file had to be 
plotted to identify its value and position. The procedure 
finally adogted was the use of a file identifier using the 
first x and y. The x was given a number like 1.03 to mark it 


as the third file of map cne and the y value was set equal 
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to the contour value. The contour value was entered as a 
large negative number, e.g. -1000, when the operator was 


unsure of the contour. 


This procedure worked well at Nurek which was digitised 
later, and hopefully the comments may save confusion if 
someone else uses this apfroach. Two methods are possible to 


reduce this digitised data into a regular grid. 


The weighted distance method is a commonly used 
technique in determining the value at a given point. The 
data suite is searched for the six, or sometimes eight or 
more, nearest neighbours to the desired point with each 
contributing to the value depending on how far they are 


away. 


There are three disadvantages to this method. The first 
is that the fall-off distance-weighting functions have to be 
chosen and these will vary with the application. Potential 
functions fall off as 1/r2 with distance, while others fall 
as 1/r3 or 1/r, but often the choice is quite arbitrary. The 
second disadvantage is that it is possible for the near 
neighbours, say six, ali to lie on one side of the required 
point, and the value will be biased toward those values 
while a seventh neighbour on the other side contributes 


nothing. This pcor choice of neighbour distribution can lead 
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to misleading results. The greatest disadvantage to the 
‘near neighbour' method lies in the cost of determining the 
neighbours. In truly random data each search point requires 
all the data points to be scanned. The Oroville data was 
fitted to a 256x256 regular grid and needed 65536 near 
neighbour searches of the 16000 data points. This method was 
considered unacceptable based on an estimate of the computer 


time necessary to extract the regular grid. 


5-1.4 Irregular Triangular Interpolation 


See eee ——— ——S ———— ae ssasac See. So eee See ee 


An alternative to the weighted distance method is to 
use triangular elements (Gold et.al. 1977). The random data 
is connected to form triangles, each being as near to 
equilateral as possible and with no overlaps. A universal 
triangle was used to enclose the data and specify the values 
at a large distance from the data. A typicai triangular 


network is shown in figure 5.1la. 


The choice of the digitising interval caused some 
difficulties since the resolution at any point was limited 
by the sampling rate. Spurs and valleys were poorly 
interpolated until extra contours were inserted at crests 
and valleys. At Oroville the data separation also had to be 


reduced when the contours were close together. 


There are two advantages to this method over the 


weighted distance method. The first is that interpolation is 
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figure 5.1: 


A typical triangular net that would be set 
up by the interfolation programme. The co- 
ordinates used are shown in (b). The three 
co-ordinates will all be positive and 
between zero and one only if the point is 
internal to the triangles. 
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always done by three points spaced around the desired point. 
The number of weighting pcints could be increased as it is 
possible to know which triangles are adjacent to a 
particular triangle since they share two mutual points. 
These neighbcuring triangles are used if continuity of the 
derivative is desired as in contouring applications. The 
main advantage to triangulation is that searching for points 
is ordered and a directed search is used to determine the 
triangle in which a particular point lies. This is made 
possible by triangular coordinates which indicate by their 
Magnitude on which side of a triangle the point lies. A 
negative coordinate, see figure 5.1, for a certain point 
indicates that the point is behind a side. The adjacent 
triangle to this side is then used to determine new 
coordinates. Only when all three triangular coordinates are 
pcsitive and between zero and one is a certain point within 
that triangle. The coordinates are then used as weighting 


functions. 


The interpolation is linear along each side and is the 
way one would probably do the problem manually - a line 
would be drawn joining each of the points and linear 
interpolation used along it. If a, b, and c are the 
triangular coordinates for a point, pf, and A, B and C the 
value at the triangle corners the weighted valued at p is: 

Z(p) = aA + bB + cC 


Gold (1977) estimates the triangular interpolation is at 
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least three times faster than the weighted distance method 
and has the added advantage of ensuring a good interpolation 


distribution. 


In the data used it was found that a lot of digitising 
could be omitted where the ccntours were equally spaced and 
parallel. Eventually about a third of the data was removed 
without any icss of accuracy. Some extra contours at crests 
and valleys were added to increase the detail in certain 


regions. 


As a guide to the cost of the method statistics were 
kept for each subset cf the Oroville data. In the case of 
set 4 (right section of the map) there were 2550 points of 
digitised data, with a further 748 points added by the 
programme tc avcid pocr interpolation when contours are 
close together. This was done by bisecting data pcints along 
the contour until the difference between points was less 


than the distance between contours. 


Searches were made fcr the values of 5358 regularly 
spaced points in this triangular net. The triangular 
arrangement requires optimisation for the triangular net so 
that better interpolaticn is possible. For this reason the 
data was inserted in sections (actually digitised files) and 
optimised after each insert. The results and costs are 


summarised in table 5.1. 


The data cbtained agreed when checked manually and the 


bos havsge ylieaps waey | eaeouse # 
bave avs aid io bates pinen 

arom. 2. +b & smote ie har 
ezeai> to 2ivesaon asa gape xoeang9 40 | anol rae 3 
nis?ae> 01 [aetsh att seeeoone sll bebbs 19% eyoliav 


sie ecitveitndd boiseu oat $b doo pdt oF sbiog Aras area } 
yo sae Sat HT ,stst Oi Lewead ae ico soadne dose 103 (sqod 
to stood 0225 stay sxodt (yee eid fo aalfooe tdpta) © tee 
add qd bebbs etutoy SPT eelpade bo omieke) «sted bomesing® 
ein g7nedooo weet uoitekogtetai te0eq Siove 92 seen 79079 | 
pHois addlog st5h tatrseced yd enob rev 22a” .sedsepot eaolo 
dust eoy | etn ioe aseutod engseeSaee eds fijnw sodssdd> odd 
seuwoseon abevsed soastetb ods pods” 
>. 2 
yiaeieper G22 to couluy de ak she az0¥ nod oz 692 _ 
seivowsizs eff .fu5 etic anit eiad ek staliog boosga. 
o@ $00 xalugabixy edt x08 nobsousedsqo aeaiupes Jnomepuss3s ui 
odo weese. abi 104 .vldtaveg 24 apitafogaetel setsed) teddy 
tine (ani2l boesskpLd yicwntyey eaoksana Ah boras0nt-seylateh 
o%6 27405 ‘bas adivass odd \Vtauepe does acsia f i 


TABLE 5.1 


206 


TRIANGULATION COST AND STATISTICS 


File Now -of Insert No. of 
Number data time optimisa- 
points (cpu) tions 
a 1289 OAS) 32 
2 577 a lieu 8 
3 Lie LOPS 4 
4 5 .018 7 
5 wS -005 ) 
6 258 2051 i 
7 22 gl OME 4 
8 80 oQES 5) 
250 ~574 70 


Number of points computed, inserted and 
optimised by the programme 


748 


S296 


|| 


TOTALS 
Searching for .5358 regular grid.points 
Time to write results 
Entering data defining the problem 


Time to read, insert and optimise 
data triangulation 


Meiers Time for 
NOt "OL insertion 
triangle & optimi- 
switches sation 


(cpu) 
6216 2.106 
2199 .826 
288 PIgd 
316 2229 
23 .079 
Tite 404 
pate 1197 
Ty 5503 
10155 4.246 
2450 0.940 
12605 5.186 
1Os027 
ky) 
AO 2 
5.186 


(cpu) 15.437 


At deferred priority on an AMDAHL-470 this 


reflects a cost of approximately $4.50. 
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technique was considered a great improvement over the manual 
method. The regular grids of Nurek and Oroville were’ used to 
compute the stresses and deflections at several times using 
the filling history of figure 3.6. The results are shown in 


figures 3.12 and 3.13. 


2-2 Discussion of Resuits 


a a ee ee ee Se ee — ——— 


The objective of this work was to see how the fiaid 
filled elastic material behaved under the influence of a 
surface load which varied both in shape and time. It is 


apparent from the figures presented that this can be done 


efficiently at any depth and time and with any loading 
history. The algebra can be solved iit three dimensions for 


Simple homogenecus material. 


There has been a large amount of numeric information 
produced and resented, but to be mneaningful it must be 
summarised into a more concise form. For this reason the 
figure 3.9, showing stresses and pressure as functions of 
time at several depths, will be analysed in detail since it 
contains all the relevant informaticn. The pressure and 
displacements fcliow curves predictable from the assumptions 


of the model. 


The displacement beccmes larger as the load increases 
and continues to increase, even though the load is kept 


constant, until it reaches a final value at a iater time. In 
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this model the water saturates the rock and, as the fluid 
flows away from the compressed zone, consolidation is 
possible even with no load changes. The initial settling in 
response to the load is a result of the elastic compression 
of the rock fabric. During the unloading there is a rebound 
Since the load is less. Gancomidimed is observed during the 
unloaded pericd. ReapEelying the lead increases the 
displacement with consolidation continuing beyond the 
refilling time. The rate of settling decreases noticeably 


several years after the refill. 


The pressure distribution is quite surprising since 
‘lobes' are observed at depth below the lake. This is 
noticed both in two and three dimensions. These regions of 
increased pressure are due to the rock and pores being 
compressed, and hence, the pressure in the pores is 
increased. The fluid cannot immediately escape from the 
pores since the motion is limited by the diffusion constant. 
The lateral extent of the lobes is governed by the fall-off 
of stress with the distance from the scurce. The top of the 
lobe depends cn the pressure boundary condition. The 
coupling choice of 0.25 causes smoother fall-off than the 
zero surface pressure condition often used in engineering 
applications which makes the lobes even more prominent (see 
figure 3.1). If the loading is Heaviside, the lobes decrease 
in amplitude as the fluid diffuses away, however, when the 


jake depth is being varied, the behaviour is more complex. 
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This is seen in figures 3.8 and 3.9 in which the fressure is 
plotted below the centre of the lake at different depths and 
times. For the lcading history chosen (figure 3.6) the 
pressure is seen to decrease while the water depth is held 
constant as expected by the diffusion process, but during 
loading the pressure depends on how fast the load is applied 
and how rapidly the fressure can diffuse away. This 
indicates how dependent the calculations are on the rate of 
loading and the rock parameters. Without an adequate 
knowledge of each at a particular dam site there is little 
value in using caiculations of this type for earthquake 


prediction. 


The centre of the Mohr circle and the maximum shear 
stress at several times after filling started are shown in 
figure 5.2 . Notice that the maximum shear stress occurs 
between 6 and 12 km deep and has a magnitude of over 1.0 bar 
except during the unload cycle. The radius of the Moby 
circle is a measure of the instability of the region. fThere 
is a close similarity in the figure 5.2 and the distribution 
of aftershocks determined at Koyna shown in figure 1.9. This 
suggests the model may provide a guide to the position of 
unstable regions, however the locations of aftershocks at 


other lakes must be known before this can be confirmed. 


The Mchr circle centre should be considered as 
important as the radius, since if two circles have the same 


shear stress, the centre position decides instability. For 
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Figure 5.2: 


Radius and centre of the Mohr circle for 
the ioading of the lake shown in figure 
3.4. The values are plotted as functions 
of depth at various stages of loading and 
unloading. The shear stresses are a 
maximum between 6. and 12 km deep directly 
below the lake. 
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this reason the locus of the point on the Mohr circle at 
which a line inclined at 30° is tangent, is plotted to see 
relative instability. The Mohr circle fcr any particular 
point can be reconstructed by projecting a line back at 60° 
to the positive stress axis. In these diagrams compressive 
stresses are positive so the figures are consistent with 
figure 1.11. Figure 5.3a indicates how the locus at 
different times was traced and 5.3b shows what is meant by 


the measure of instability for a certain Mohr circle. 


Since tectonic stresses are present the origins of ail 
these diagrams represents the tangent point for the tectcnic 
stresses, and motions cf the total Mohr circle are all 
relative to this point. The larger the instability, the 
larger is the deviation from t:she initial stressed condition, 
and supposedly, the area is more likely to fail. it is 
necessary to assume that the area is highly prestressed and 
close to failure with the anomalous stresses acting as 
triggers. The locus diagrams will differ depending on the 


assumed stress envircnment and faulting mechanism expected. 


The loci of four points were plctted at each of the 
depths of 3, 6 and 12 km. Belcw 12 km most effects decrease 
and need not be shown since the stresses are always larger 
elsewhere. Fcur fositions ‘'1', Lo, *3* and ‘'4* were 
examined at each of the depths with position 1 directly 
under the lake and ‘4* off to the side. These positions are 


shown in the figure 5.4 . The lcci of the 12 selected points 
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Figure 5.3: 


The method of determining the locus of the 
Mohr circle with increasing time is shown 
in (a). The “instability" measure is shown 
as a projection of the Mohr circle onto a 
line at 60° centered on the initial 
prestressed condition. The larger the 
instability the larger the triggering 
force available at that time and place. 
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Figure 5.4: Positions fcr which the loci of the Mohr 
circle were plotted in figure Seite 
Position 1 is directly under the lake and 
4 is offset by about 40 km from. the 
centre. 
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are shown in figure 5.5 and several of the paths are very 
unusual. These loci are for the normal faulting stress 
convention. The detail is felt to be real since the curves 
were generated using two time intervals (0.1-4.yrs and 4.- 
6.yrs), each sampled at 51 log-spaced points. A zero radius 
indicates that the X and Z components of the stress are 
equal and movement to the left and right is mostly a result 
of fluid pressure changes as illustrated by the simple 
models of figure 1.11. Since the pressure increases during 
compressions resulting frcm the loading, then diffusion will 
always move the circle towards the stabilising direction. 
This does nct however prevent the radius of the Mohr circle 
becoming larger as time prcegresses, but the general result 
observed is that after the refilling perturbation, the rock 
is always more stable. The figure 5.5 attempts to summarize 
most of the results of the research since it shows the 


seismic risk at particular regions for various times. 


At 1 year the most unstable region is at six kilometers 
depth immediately below the lake. The instability is about 
1.2 bars aktove the tectonic level and, depending on the 


sensitivity of the area, may trigger an earthquake. 


At 2 years, at the end of the first filling, the most 
unstable position is at 12 km directly below the lake, with 
all other regions having a ievel of instability less than 
this. This position at two years is more unstable than it 


waS at one year. The stresses at two years, at depth 12 kam 
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Figure 5.5: 


Loci of the Mohr circle for the points 
indicated in figure 5.4. The closer a 
point is tc a line y= tan (30°) *x+C the 
more unstable is the region during the 
initial filling At some positions 
instability is at 1 while at others it is 
at 2 years. Often the largest instability 
occurs at 5.5 (after the refill). The 
water level followed the history shown in 
figure 3.4. 
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below the jake are just slightly larger than those 
experienced earlier at 6 km one year after filling stated. 
If an event hadn't occurred at one year and 6 km deep then 
it may happen at 12 km and 2 years. The stresses at the 
positions 2, 3, and 4 are much smaller than directly under 
the lake and failure most probably will be initiated 
directly under the lake at position 1 where the instability 


is much larger. 


Between two and four years while the lake's depth is 
being kept constant the region stabilises as the high 


Fressures diffuse away. 


It is possible to conclude that the initial filling is 
associated with an unstakle period. This is consistent with 
the figure 1.10. Both Vajont and Kremasta had onset of 
seismicity during the first loading. Unfortunately the exact 
positions of epicentes are unknown and it is impossible to 
check if they lay directly under the deepest section of the 


lakes. 


The next period of the lake history involves an 
emptying, followed later by a refilling to 10 % deeper than 
earlier. This type of history was included to see what 
connection it has, if any, with the fossible delayed 
earthquake at Oroville which occurred just after a refiil 


such as this. 


At the 3 and 6 km depths below the lake the anomalous 
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stresses at 5.5 years are less than those experienced at an 
earlier time. It must be concluded that these areas would 
not have earthquakes’ since they were stressed to a higher 
level at an earlier time. At 12 km below the lake an 
interesting result is observed. The stresses after refilling 
are larger than at any previous time or position. If no 
earthquakes had occurred to relieve stresses, the refilling 
is the time at which earthquakes are most likely to happen. 
After unloads and reloads a marked increase in the 
seismicity is observed, always starting after the water 
depth starts to increase. The calculations agree with this 
Observation. The timing of tne event at Oroville could 
perhaps have been predicted from this model, however, the 
off-set must be explained on geological grounds as the off- 
set is larger than predicted. The diagrams showing loci have 
been plotted with a stress sign conventicn consistant with 


the normal faulting cbserved at Oroville. 


The stresses are always largest immediately under the 
lake, but on the assumpticn that this area will have an 
earthquake which will release the locally stored strains, 
the question becomes which area to the sides are most 
unstable? To determine this, the ‘instability' has been 
plotted at different times on figure 5.6. It is apparent 
that there is a zone where the x and z effective stresses 
are equal . To either side of the central area are regions 


with negative instability - these areas are more stable than 
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Figure 5.6: 


Instability (see figure 5.4) at different 
times during the loading. The results are 
shown to 12 km depth with instabilities 
(stabilised zones) having valves less than 
zero nct being contoured. 
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before filling. 


Locus were redrawn using the stress convention in a 
thrust faulting environment with the vertical stress being 
smaller than the hcrizcntal. The same stress values were 
used to determine the locus shown in figure 5.7. It can be 
seen that the conclusions discussed abcve no longer apply. 
The initial filling leads to instability only off to the 
sides of the lake. The instability levels are very small and 
if failure did not cccur during the initial filling, then it 
would appear unlikely to happen later.From the figure 5.6 it 
appears that the instability may be large at shallow depths 


at the end of the down-draw in a thrust environment. 


To illustrate the effect of changing the fluid 
parameters , a ncrmal faulting locus plot was made using 
three time intervals (0.01-4yrs, 4-6yrs,6.-25yrs) fea 
alpha=1 and Q cf 5.Mbar. This is shown in figure 5.8. Again 
the largest instability is at five years or later for the 
locus found. The movement of the circle to the left 
indicates that the pressure is increasing with time so the 
effective stresses are decreasing. The results at 25 years 
are very close to the final locus position. This illustrates 
that the conclusions can vary considerably if the fluid 


parameters are uncertain. 


3 Concluding Remarks 
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Figure 5.8: 


Instability in a normal faulting area with 
the Q reduced so that the diffusion may 
take place. The coupling is one in this 
case. The result differs considerably from 
those shown earlier. 
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it appears that the model can explain many of the 
observed seismic effects at reservoirs. The stresses, pore 
pressures and displacement functions each have terms 
depending on the convolution of the rate of change of water 
level, on the water depth, and on a diffusion time. The 
instability during initial loading and again during an 


unload-reload situation is encouraging. 


Applicaticn of this werk to a particular reservoir in 
detail has been avoided for three reasons. Firstly, not 
enough is krown about the geologic and elastic behaviour of 
the rocks; secondly, nothing is  kncwn about the initial 
stresses and hydrology; and finally, the epicentral 
positions and depths are usually inaccurately known. This 
leads tc the question of future research. How can the 
peaultt= and conclusions made here be tested? Several 
experiments can be made to check the results, however, each 


is expensive and are major projects. 


Of utmost importance is a test to find the amount of 
triggering stresses needed to open fractures in the 
reservoir area and this cculd be done by hydrofracturing at 
a deep test hole near the site. The test hole could also 
provide geological samples, insitu permeability results and 
information on the water table, both top and bottom. A 
seismic velocity and resistivity log down the hole would 
also be useful. The resistivity log could be used tc monitor 


Ghanges in the water table cheaply and easily during and 
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after filling. Several shalicw holes could be drilled around 
the test site, again tc mcnitor water table changes, but 
also as source holes in diffusion studies over long 


distances so that a 'bulk permeability" is obtained. 


It is important that the hydrology and fluid parameters 
be known. This is not a simple task over the large scales 
discussed here. Regicnal hydrology would have to be obtained 
if initial conditions are to be imposed on the fluid 
behaviour of the rock. A guide to the rock parameters may be 
obtained by modelling the behaviour of eee reservoirs in 


the region. 


A promising line of predictive research lies in the 
monitoring of the delay time of a seismic Signal over a 
fixed baseline. If sufficient acoustic energy can be pumped 
into the ground, it may be possible to pass the waves 
beneath the reservoir at several kilometers depth. Since the 
Material properties are stress dependent, an increase of the 
stress is thought by scme researchers (Reasonberg and Aki, 
1974) to change the velocity sufficiently so that one bar 
stress changes may be detected with good seismic timing. The 
amount of delay is uncertain, but if it is sufficient, 
insitu stresses could be continuously monitored. The seismic 
energy could perhaps be obtained from construction blasts. 
Strain meters may also be useful in monitoring strain 
accumulation. Eventually an extensive study will have to be 


done into the depths and lateral migrations of earthquakes 
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at dams so that models of the type presented here can be 


compared. Hopefully, scme of the data being collected now 
will be soon available. The research groups at Nurek and 
Lamont are currently cbhtaining good seismic locations. In 
any event, an array of seismometers should be placed at all 


new dam sites to monitor any activity that occurs. 


Although the mathematical calculations done here are 
gratifyingly close to what little observational data exists, 
the reader is cautioned that application of mathematical 
models to real problems in earth science is more of an 
uncertain art than an unambiguous science. Mathematical 
modelling like this work indicates trends, and provides 
Gguantitative measures which can make the weight of evidence 
for, or against, seismic risk at a dam site far more 
convincing. it is not by itself sufficient evidence, but it 


is an essential compcnent of the evidence. 
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